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I.  Introduction 

This  Postdoctoral  Training  Award  (W81XWH-08-1-0127,  entitled  “Accurate  and  Fast 
Localization  of  Prostate  for  External  Beam  Radiation  Therapy”)  was  awarded  to  the  principal 
investigator  (PI)  for  the  period  of  February  15,  2008  —  February  14,  2010.  This  is  the  annual 
report  for  the  first  funding  period  (February  15,  2008  —  February  14,  2009).  The  purpose  of  this 
project  is  to  develop  a  tomosynthesis-based  method  for  fast  and  accurate  prostate  localization 
during  IMRT.  The  specific  aims  of  this  project  are:  1)  To  determine  the  protocol  for  acquiring 
projection  data  for  tomosynthesis,  2)  to  determine  the  optimal  tomosynthetic  image 
reconstruction  algorithm,  3)  to  determine  an  accurate  non-rigid  registration  method  for 
registering  the  daily  tomosynthetic  images  to  the  planning  CT  images  of  the  prostate,  4)  to 
develop  and  evaluate  the  tomosynthesis-based  repositioning  protocols.  Under  the  generous 
support  from  the  U.S.  Army  Medical  Research  and  Materiel  Command  (USAMRMC),  the  PI  has 
contributed  significantly  to  the  field  of  prostate  cancer  research  by  applying  the  physics  and 
engineering  knowledge.  A  number  of  conference  abstracts  and  refereed  journal  publications  have 
been  resulted  from  the  support.  In  this  report,  the  past  year’s  research  activities  of  the  PI  are 
highlighted. 

II.  Body 

During  the  past  year,  the  PI  has  successfully  implemented  a  ray-tracing  method  to  generate  the 
projection  image  from  a  digital  mathematical  phantom  and  used  Trilogy  oncology  system  to 
acquire  projection  data  of  an  anthropomorphic  phantom  with  various  protocols.  The  PI  has  also 
implemented  tomographic  image  reconstruction  algorithms  including  both  analytical  filtered 
back-projection  (FBP)-type  method  and  statistics-based  iterative  image  reconstruction  algorithm. 
The  analytical  image  reconstruction  algorithm  is  based  on  the  widely-used  Feldkamp  Davis  and 
Kress  (FDK)  algorithm.  To  enhance  the  image  of  reconstructed  tomographic  image,  a  statistics- 
base  sinogram  smoothing  algorithm  has  been  developed.  Two  journal  papers  (ref.l  and  2)  and  a 
conference  proceeding  (ref  3)  have  been  published  based  on  the  developed  algorithm.  The 
iterative  image  reconstruction  algorithm  is  based  on  statistical  properties  of  measured  projection 
data,  where  the  noise  modeling  of  measured  data  is  incorporated  into  the  penalized  weighted 
least-squares  objective  function.  To  improve  the  image  resolution  of  reconstructed  image,  a  new 
prior  model  is  proposed.  One  conference  abstract  (ref.4)  and  one  journal  paper  (ref  5)  have  been 
published  based  the  proposed  algorithm.  Currently,  the  PI  is  conducting  research  on  image 
reconstruction  and  registration  for  tomosynthesis  using  few  projections  from  the  cone-beam  CT. 
One  paper  regarding  the  number  of  projection  required  for  accurate  tomosynthetic  image 
reconstruction  and  registration  is  under  preparation. 

III.  Key  Research  Accomplishments 

•  Developed  a  simulation  package:  including  projector  to  generate  projection  image  from  a 
digital  mathematical  phantom  and  FDK  reconstruction  algorithm  to  reconstruct 
tomosynthesis  image. 

•  Developed  a  statistics  -based  sinogram  smoothing  algorithm  to  suppress  noise  in 
projection  data  and  improve  image  quality  in  reconstructed  image. 
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•  Developed  an  iterative  image  reconstruction  algorithm  with  a  new  prior  constraint  to 
improve  image  quality  and  resolution  in  tomographic  images. 

IV.  Reportable  Outcomes 

The  following  is  a  list  of  publications  resulted  from  the  grant  support  in  the  last  funding  period. 

Refereed  Publications: 

1.  J.  Wang,  L.  Zhu,  and  L.  Xing,  “Noise  reduction  in  low-dose  X-ray  fluoroscopy  for 
Image  Guided  Radiation  Therapy  (IGRT)”,  Int  J  Radiat  Oncol  Biol  Phys,  in  press, 
2009 

2.  J.  Wang,  T.  Li,  and  L.  Xing,  “Iterative  image  reconstruction  for  CBCT  using  edge¬ 
preserving  prior”.  Medical  Physics,  vol.  36,  pp.  252-260,  2009 

3.  L.  Zhu,  J.  Wang,  and  L.  Xing,  “Noise  suppression  in  scatter  correction  for  Cone- 
Beam  CT”,  Medical  Physics,  vol.  36,  pp.  741-752,  2009 

4.  J,  Wang,  T.  Li,  Z.  Liang  and  L.  Xing,  “Dose  reduction  for  kilovoltage  cone-beam 
computed  tomography  in  radiation  therapy”.  Physics  in  Medicine  and  Biology,  vol. 
53,  pp.  2897-2909,  2008 

Conference  Abstract: 

1.  J.  Wang,  A.  Chai,  L.  Xing,  “Noise  correlation  in  CBCT  projection  data  and  its 
application  for  noise  reduction  in  low-dose  CBCT”,  poster  presentation  in  2009 
SPIE  Medical  Imaging  conference,  Orlando,  FL 

2.  J.  Wang,  T.  Li,  and  L.  Xing,  “Low-dose  CBCT  Imaging  for  External  Beam 
Radiotherapy”,  oral  presentation  in  2008  ASTRO  Annual  Meeting,  Boston,  MA 

3.  X.  Zhang,  J.  Wang,  L.  Zhu,  and  L.  Xing,  “Low-dose  X-ray  fluoroscopy  for  Image 
Guided  Radiation  Therapy  (IGRT)”,  poster  presentation  in  2008  ASTRO  Annual 
Meeting,  Boston,  MA 

4.  J.  Wang,  T.  Li,  Z.  Liang  and  L.  Xing,  “Dose  reduction  for  kilovoltage  cone-beam 
computed  tomography  in  radiation  therapy”,  oral  presentation  in  2008  AAPM 
Annual  Meeting,  (selected  for  long  presentation  at  the  John  S,  Laughlin  Science 
Council  Research  Symposium),  Houston,  TX 

5.  J.  Wang,  L.  Zhu,  A.  Chai,  and  L.  Xing,  “Temporal  filtering  of  noise  in  low-dose  x- 
ray  fluoroscopy”,  poster  presentation  in  2008  AAPM  Annual  Meeting,  Houston,  TX 

6.  L.  Zhu,  J,  Wang,  Y.  Xie,  J.  Starman,  R.  Fahrig,  and  L.  Xing,  “A  patient  set-up 
protocol  based  on  partially  blocked  cone-beam  CT”,  poster  presentation  in  2008 
AAPM  Annual  Meeting,  Houston,  TX 

7.  J.  Wang,  T.  Li,  and  L.  Xing,  “Iterative  image  reconstruction  for  on-board  CBCT”, 
poster  presentation  in  2008  Electronic  Portal  Imaging  &  Positioning  Devices,  San 
Francisco,  CA 

V.  Conclusions 

In  summary,  an  infrastructure  has  been  established  to  execute  the  proposed  research.  Novel 

image  reconstruction  and  processing  algorithms  have  been  proposed  for  the  treatment  of  prostate 

cancer.  A  few  milestones  have  been  achieved  toward  the  general  goal  of  the  project. 

Implementation  and  evaluation  of  the  image  registration  algorithms  are  current  underway  to 
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integrate  the  developed  image  reconstruction  algorithms  to  improve  the  current  prostate  IMRT 
treatment. 

VI.  Reference: 

1 .  L.  Zhu,  J.  Wang,  and  L.  Xing,  “Noise  suppression  in  scatter  correction  for  Cone- 
Beam  CT”,  Medical  Physics,  vol.  36,  pp.  741-752,  2009 

2.  J.  Wang,  T.  Li,  Z.  Liang  and  L.  Xing,  “Dose  reduction  for  kilovoltage  cone -beam 
computed  tomography  in  radiation  therapy”.  Physics  in  Medicine  and  Bioloyy, 
vol.  53,  pp.  2897-2909,  2008 

3.  J.  Wang,  A.  Chai,  L.  Xing,  “Noise  correlation  in  CBCT  projection  data  and  its 
application  for  noise  reduction  in  low-dose  CBCT”,  poster  presentation  in  2009 
SPIE  Medical  Imaging  conference,  Orlando,  FL 

4.  J.  Wang,  T.  Li,  and  L.  Xing,  “Low-dose  CBCT  Imaging  for  External  Beam 
Radiotherapy”,  oral  presentation  in  2008  ASTRO  Annual  Meeting,  Boston,  MA 

5.  J.  Wang,  T.  Li,  and  L.  Xing,  “Iterative  image  reconstruction  for  CBCT  using 
edge-preserving  prior”.  Medical  Physics,  vol.  36,  pp.  252-260,  2009 

VII.  Appendices 

Refereed  publications  in  the  last  funding  period 
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Abstract 

Kilovotage  cone-beam  computed  tomography  (kV-CBCT)  has  shown  potentials 
to  improve  the  accuracy  of  a  patient  setup  in  radiotherapy.  However,  daily  and 
repeated  use  of  CBCT  will  deliver  high  extra  radiation  doses  to  patients.  One 
way  to  reduce  the  patient  dose  is  to  lower  mAs  when  acquiring  projection  data. 
This,  however,  degrades  the  quality  of  low  mAs  CBCT  images  dramatically 
due  to  excessive  noises.  In  this  work,  we  aim  to  improve  the  CBCT  image 
quality  from  low  mAs  scans.  Based  on  the  measured  noise  properties  of  the 
sinogram,  a  penalized  weighted  least-squares  (PWLS)  objective  function  was 
constructed,  and  the  ideal  sinogram  was  then  estimated  by  minimizing  the 
PWLS  objection  function.  To  preserve  edge  information  in  the  projection 
data,  an  anisotropic  penalty  term  was  designed  using  the  intensity  difference 
between  neighboring  pixels.  The  effectiveness  of  the  presented  algorithm  was 
demonstrated  by  two  experimental  phantom  studies.  Noise  in  the  reconstructed 
CBCT  image  acquired  with  a  low  mAs  protocol  was  greatly  suppressed  after 
the  proposed  sinogram  domain  image  processing,  without  noticeable  sacrifice 
of  the  spatial  resolution. 

(Some  figures  in  this  article  are  in  colour  only  in  the  electronic  version) 


1.  Introduction 

Integration  of  the  kilovotage  cone-beam  computed  tomography  (kV-CBCT)  with  a  linear 
accelerator  makes  it  possible  to  acquire  a  high-resolution  volumetric  image  of  a  patient  at 
a  treatment  position.  There  is  growing  interest  in  using  on-board  kV-CBCT  for  a  patient 
treatment  position  setup  and  dose  reconstruction  in  radiotherapy  (Xing  et  al  2006,  Yang  et  al 
2007,  Lee  et  al  2008).  However,  the  repeated  use  of  kV-CBCT  during  the  course  of  a  treatment 
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has  raised  concerns  of  an  extra  radiation  dose  delivered  to  patients  (Brenner  and  Hall  2007, 
Islam  et  al  2006,  Wen  et  al  2007).  It  has  been  reported  (Wen  et  al  2007)  that  the  dose  delivered 
from  Varian’s  kV-CBCT  system  with  current  clinical  protocols  is  more  than  3  cGy  for  central 
tissue  and  about  5  cGy  for  most  of  the  peripheral  tissues  during  an  IMRT  (intensity-modulated 
radiation  therapy)  treatment  course  for  prostate  cancer.  The  extra  radiation  exposure  to  normal 
tissue  during  kV-CBCT  will  significantly  increase  the  probability  of  stochastic  risk  of  inducing 
cancer  and  genetic  defects.  Based  on  the  ALARA  (as  low  as  reasonably  achievable)  principle, 
the  unwanted  kV-CBCT  dose  should  be  minimized  to  fully  realize  its  advantages  of  precise 
target  localization  during  radiotherapy  (Murphy  et  al  2007). 

One  way  to  reduce  the  radiation  dose  delivered  to  patients  during  the  kV-CBCT  procedure 
is  to  acquire  CT  projection  data  with  a  lower  mAs  level  (can  be  realized  by  reducing  the 
tube  current  or  pulse  time).  However,  the  image  quality  of  the  projection  image  and  the 
reconstructed  CBCT  image  will  be  degraded  due  to  excessive  quantum  noise  as  a  result  of 
a  low  mAs  protocol.  Conventionally,  noise  in  CT  is  suppressed  by  using  a  low-pass  filter 
to  attenuate  the  high-frequency  component  of  the  projection  data  during  reconstruction.  The 
high-frequency  component  contains  information  of  both  noise  and  edges,  where  a  simple 
low-pass  filter  cannot  differentiate  edge  information  from  noise.  Therefore,  noise  reduction 
using  a  low-pass  filter  will  result  in  loss  of  edges,  which  is  not  desirable  for  CT  imaging. 
Several  edge-preserving  filters  (Hsieh  1998,  Kachelriess  et  al  2001,  Zhong  et  al  2004)  have 
been  proposed  to  reduce  noise  in  CT  images  based  on  local  characteristics  of  the  projection 
data  elements.  More  recently,  statistics-based  image  domain  (Li  et  al  2005a)  and  sinogram 
domain  restoration  algorithms  (Li  et  al  2004,  La  Riviere  2005,  La  Riviere  and  Billmire  2005, 
Wang  et  al  2006)  have  shown  advantages  in  noise  reduction  and  edge  preservation  for  low- 
dose  fan-beam  CT.  In  the  meantime,  noise  properties  of  CT  projection  data  have  been  under 
investigation  (Li  et  al  2004,  Whiting  et  al  2006)  and  the  noise  model  of  the  sinogram  data  in 
Radon  space  (i.e.  line  integrals)  has  been  validated  by  experimental  studies  (Wang  et  al  2008). 
In  this  work,  we  aim  to  improve  the  low-dose  CBCT  image  quality  by  reducing  noise  in  the 
CBCT  sinogram  before  image  reconstruction.  The  noise  reduction  algorithm  incorporates 
the  noise  modeling  of  the  CT  sinogram  data  in  Radon  space  (line  integrals)  to  construct 
a  penalized  weighted  least-squares  (PWLS)  objective  function  (Fessler  1994,  Sukovic  and 
Clinthorne  2000).  The  ideal  solution  of  the  line  integrals  is  then  estimated  by  minimizing 
the  PWLS  objective  function.  The  weighted  least  square  is  based  on  the  first  and  second 
moments  of  the  noise  in  the  sinogram  data  and  an  anisotropic  penalty  is  designed  to  preserve 
the  edges  in  the  sinogram.  CBCT  images  are  reconstructed  by  using  the  Feldkamp-Davis- 
Kress  (FDK)  (Feldkamp  et  al  1984)  algorithm  after  all  sinogram  images  are  processed  by  the 
PWLS  criterion  sequentially.  The  effectiveness  of  the  PWLS-based  noise  reduction  algorithm 
is  demonstrated  by  two  experimental  phantom  studies. 


2.  Methods  and  materials 


2.7.  CBCT  sinogram  smoothing 


Ideally,  the  line  integral  of  attenuation  coefficients  can  be  calculated  by 


=  In 


Ni 


(1) 


where  Wo  and  W,  are  the  incident  photon  number  and  detected  photon  number  at  the  detector 
bin  i  respectively.  For  ease  of  presentation,  we  refer  the  measurement  as  a  photon  number.  In  a 
real  x-ray  CBCT  system,  the  measured  signal  is  total  energy  deposit  on  the  flat -panel  detector. 
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In  the  following  of  this  paper,  we  refer  the  value  of  p,  as  the  sinogram  datum  at  the  detector 
bin  i.  Mathematically,  the  PWLS  cost  function  in  the  sinogram  domain  can  be  written  as 

^{p)^{y-P)'^^-\y-P)  +  pR(p).  (2) 

The  hrst  term  in  equation  (2)  is  a  weighted  least-squares  criterion,  where  y  is  the  vector  of 
the  measured  sinogram  data  and  p  is  the  vector  of  the  ideal  sinogram  data  to  be  estimated. 
The  symbol  T  denotes  the  transpose  operator.  The  matrix  E  is  a  diagonal  matrix  and  its  /th 
element  is  the  variance  of  sinogram  data  at  the  detector  bin  i.  The  second  term  in  equation 
(2)  is  a  smoothness  penalty  or  a  priori  constraint,  where  /J  is  the  smoothing  parameter  which 
controls  the  degree  of  agreement  between  the  estimated  and  the  measured  data. 

The  element  of  the  diagonal  matrix  E  is  the  variance  of  the  corresponding  sinogram 
datum,  and  it  determines  the  contribution  of  each  sinogram  datum  to  the  cost  function.  Based 
on  the  sinogram  noise  modeling  in  Li  et  al  (2004)  and  Wang  et  al  (2008),  the  variance  of  the 
sinogram  datum  can  be  estimated  by 

af  =  exp(p,)/Wo-  (3) 

For  a  hxed  incident  photon  number  A^,o,  a  sinogram  datum  with  a  larger  value  will  have  a 
larger  variance  and  therefore  less  contribution  to  the  cost  function  since  the  weight  of  that 
measured  datum  is  l/ af  as  dehned  in  equation  (2).  This  can  be  understood  by  the  following 
observation.  A  larger  sinogram  datum  value  p,  at  the  detector  bin  i  indicates  less  x-ray 
photons  being  detected,  i.e.  smaller  A,  in  equation  (1),  or  more  photons  being  attenuated 
along  the  projection  path  i.  A  detector  bin  receiving  less  photons  will  be  associated  with  a 
smaller  signal-to-noise  ratio  (SNR)  based  on  the  Poisson  noise  nature  of  the  detected  x-ray 
photons.  Therefore,  the  weighted  least-squares  criterion  reflects  the  above  observation  that 
the  measured  datum  with  a  lower  SNR  will  contribute  less  for  estimation  of  its  ideal  sinogram 
datum. 

To  calculate  the  sinogram  datum  variance  at  the  detector  bin  i  via  equation  (3),  we  need 
to  estimate  the  incident  photon  number  A,o  for  calculation  of  the  sinogram  variance.  The 
incident  photon  number  is  mainly  determined  by  the  protocols  of  tube  current  and  the  duration 
of  x-ray  pulse  (i.e.  mAs).  Ideally,  the  incident  x-ray  flux  from  the  tube  would  be  calibrated 
as  uniform  as  possible  across  a  held  of  view  (FOV),  i.e.  A,o  is  a  constant  for  all  the  detector 
bins.  In  reality,  the  x-ray  flux  is  modulated  to  consider  the  concavity  shape  of  the  human 
body  by  the  bow-tie  attenuating  Alter  prior  to  arrival  at  the  patients.  Therefore,  the  incident 
photon  number  will  not  be  a  constant  across  the  FOV.  To  estimate  the  incident  intensity  over 
the  FOV  at  a  specific  mAs  level,  we  performed  the  air  scan  and  then  averaged  the  projection 
image  over  all  projection  view  angles.  Figure  1  shows  the  incident  x-ray  intensity  with  the 
tube  current  80  mA  and  duration  of  pulse  10  ms.  The  incident  x-ray  intensity  can  then  be  used 
for  estimation  of  the  sinogram  data  variance  {a^f}- 

The  penalty  term  in  equation  (2)  is  a  prior  or  smoothing  constraint,  which  encourages  the 
equivalence  between  neighboring  data  elements.  In  Li  et  al  (2004)  and  Wang  et  al  (2006),  a 
penalty  of  a  quadratic  form  with  equal  weights  for  all  neighbors  has  been  used  for  sinogram 
smoothing  of  fan-beam  CT; 

Rip)  =  ~  i"^) 

n 

where  n  represents  four  nearest  neighbors  around  pixel  i  and  is  the  weight  for  neighbor 
n.  With  an  equal  weight  for  the  four  nearest  neighbors,  these  neighbors  play  an  equivalent 
role  in  constraining  the  solution.  As  such,  it  provides  a  uniform  regularization  without 
considering  details  of  intensity  variation  and  possibly  the  presence  of  edges  in  the  sinogram 
image.  To  preserve  the  edge  information  in  the  sinogram  image  of  CBCT,  we  propose  to  use 
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Figure  1.  Incident  x-ray  intensities  across  the  field  of  view  with  80  mA  tube  current  and  10  ms 
pulse  time.  Relative  intensity  is  mainly  caused  by  the  bow-tie  filter. 


anisotropic  weights  for  different  neighbors  in  the  sinogram  image.  The  weight  of  the  neighbor 
is  determined  by  the  magnitude  of  difference  between  neighbors  and  the  concerned  pixel.  For 
a  larger  difference  between  the  neighbor  and  the  pixel,  the  coupling  between  them  should  be 
weaker  and  the  weight  w,,,  should  be  smaller.  This  form  of  weight  is  chosen  the  same  as  the 
conducting  coefficient  in  the  well-known  anisotropic  diffusion  filter  (Perona  and  Malik  1990): 


(5) 


Win  =  exp 


where  the  gradient  determines  the  strength  of  the  diffusion  during  each  iteration  and  the 
parameter  &  was  chosen  as  90%  of  histogram  of  the  gradient  magnitude  of  the  sinogram  to  be 
processed  (Perona  and  Malik  1990). 

Minimization  of  the  objective  function  2  can  be  performed  efficiently  by  the  iterative 
Gauss-Seidel  updating  strategy.  The  updating  formula  for  the  solution  of  p  is  given  by 


(6) 


Pi 


where  the  index  k  represents  the  iterative  number,  denotes  those  two  nearest  neighbors  of 
i  whose  index  is  smaller  than  /,  77?  denotes  those  two  nearest  neighbors  of  i  whose  index  is 


larger  than  i  and  A^^denotes  these  four  nearest  neighbors  of  pixel  i  in  the  sinogram  image.  The 
initial  of  p  is  given  by  the  measured  data  y. 

2.2.  On-board  kV-CBCT 

The  cone-beam  CT  projection  data  were  acquired  by  Exact  Arms  (kV  source/detector  arms) 
of  a  Trilogy(tm)  treatment  system  (Varian  Medical  Systems,  Palo  Alto,  CA).  The  number  of 
projections  for  a  full  360°  rotation  is  around  634.  The  dimension  of  each  acquired  projection 
image  is  397  mm  x  298  mm,  containing  1024  x  768  pixels.  The  system  has  a  FOV  of  25  cmx 
25  cm  (full-fan  mode)  in  the  transverse  plane  and  17  cm  in  the  longitudinal  direction,  which 
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Figure  2.  Illustration  of  the  anthropomorphic  head  phantom  used  for  evaluation  of  the  PWLS 
algorithm. 


can  be  increased  to  45  cm  x  45  cm  in  the  transverse  plane  by  shifting  the  detector  laterally 
(half-fan  mode). 

Two  phantoms  were  used  to  evaluate  the  performance  of  the  proposed  PWLS  algorithm 
in  this  study.  The  hrst  phantom  is  a  commercial  calibration  phantom  CatPhan®  600  (The 
Phantom  Laboratory,  Inc.,  Salem,  NY).  Details  about  the  CatPhan®  600  phantom  can  be 
found  in  Li  et  al  (2005a).  The  second  one  is  an  anthropomorphic  head  phantom  (see  hgure  2). 
For  each  phantom,  the  x-ray  tube  current  was  set  at  10  mA  (low  dose)  and  80  mA  (high  dose) 
during  acquisition  of  CBCT  projection  images.  At  both  mA  levels,  the  duration  of  the  x-ray 
pulse  at  each  projection  view  was  10  ms.  The  tube  voltage  was  set  to  125  kVp  during  all  data 
acquisitions.  After  each  sinogram  acquired  with  the  low-mAs  protocol  was  processed  by  the 
PWLS  algorithm  described  above,  the  CBCT  image  was  reconstructed  by  the  FDK  algorithm. 
The  voxel  size  in  the  reconstructed  image  is  0.5  x  05  x  0.5  mm^. 


3.  Results 

3.1.  CatPhan®  600  phantom 

We  first  tested  the  proposed  algorithm  on  the  CatPhan®  600  phantom.  Several  representative 
slices  of  the  reconstructed  CBCT  are  shown  in  hgures  3, 4  and  6.  In  each  of  these  hgures,  (a)  is 
the  FDK  reconstructed  image  from  the  projection  data  acquired  with  10  mA  tube  current,  (b)  is 
the  FDK  reconstructed  image  from  the  sinogram  processed  by  the  proposed  PWLS  sinogram 
smoothing  algorithm  and  (c)  is  the  FDK  reconstructed  image  from  the  sinogram  obtained  with 
80  mA  tube  current. 

Figure  3  shows  that  one  slice  of  image  contains  a  point-like  object,  which  mimics  a 
hducial  marker.  In  figure  3(a),  the  point  source  is  difficult  to  be  observed.  After  the  sinogram 
was  processed  by  the  PWLS  algorithm,  the  reconstructed  image  (figure  3(b))  is  very  similar 
to  that  obtained  with  a  high  mA  protocol  (hgure  3(c)).  The  point  source  was  well  recovered 
and  easy  to  be  detected. 
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Figure  3.  One  slice  of  the  FDK  reconstructed  image  of  the  CatPhan®  600  phantom  containing 
a  point-like  object:  (a)  from  projection  images  acquired  with  10  mA  tube  current,  (b)  after  the 
sinogram  acquired  with  10  mA  tube  cuirent  is  processed  by  the  PWLS  algorithm  and  (c)  from 
projection  images  acquired  with  80  mA  tube  current. 


Figure  4  shows  that  one  slice  of  image  contains  several  strips  with  different  sizes  and 
contrasts,  which  can  be  used  to  study  the  edge  information  in  the  reconstructed  images.  The 
CT  image  reconstructed  from  the  PWLS-processed  sinogram  is  comparable  to  that  obtained 
with  the  80  mA  protocol  in  terms  of  detectability  of  the  strips;  see  ROI2  in  figure  4(c).  To 
show  the  difference  between  figures  4(a),  (b)  and  (c),  in  figure  5  we  plotted  horizontal  profiles 
along  the  central  strips  (see  ROll  in  figure  4(c)).  It  can  be  observed  that  the  edges  are  well 
preserved  (compare  profiles  through  figures  3(b)  and  (c)),  while  noise  is  effectively  suppressed 
(compare  profiles  through  figures  3(a)  and  (b)). 

To  further  quantitatively  evaluate  the  effectiveness  of  the  PWLS  sinogram  smoothing 
algorithm,  we  calculated  the  contrast-to-noise  ratio  (CNR)  at  different  regions  of  interest 
(ROIs)  in  the  images  shown  in  figure  6.  The  CNR  is  defined  as 


CNR  = 


I  Ms  Mb  I 


(7) 
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Figure  4.  One  slice  of  the  FDK  reconstructed  image  of  the  CatPhan®  600  phantom  containing 
several  strips:  (a)  from  projection  images  acquired  with  10  mA  tube  current,  (b)  after  the  sinogram 
acquired  with  10  mA  tube  cument  is  processed  by  the  PWLS  algorithm  and  (c)  from  projection 
images  acquired  with  80  mA  tube  current. 


Table  1.  CNRs  of  four  ROIs  in  figure  5. 


ROIl 

ROE 

ROE 

ROM 

ROE 

80  mA 

1.83 

7.31 

4.75 

1.51 

0.89 

10  mA 

0.82 

2.70 

1.68 

0.49 

0.36 

PWLS  10  mA  ,6  =  0.05 

1.92 

6.88 

4.75 

1.33 

0.85 

where  /Xs  is  the  mean  value  of  the  signal  and  /Xb  is  the  mean  value  of  the  background.  Five 
circular  objects  (indicated  by  arrows  in  figure  6)  with  different  intensities  were  chosen  to 
calculate  CNRs.  Table  1  lists  the  CNRs  of  these  five  ROIs.  After  a  10  mA  sinogram  was 
processed  by  the  PWLS  algorithm,  the  CNR  in  the  reconstructed  image  improved  significantly. 
It  can  be  observed  that  the  CNR  of  a  PWLS-processed  10  mA  image  is  comparable  to  that  of 
the  image  acquired  with  the  80  mA  protocol. 
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Figure  5.  Profiles  through  the  central  strips  in  figure  4  (indicated  by  ROIl). 


3.2.  Anthropomorphic  head  phantom 

Results  of  the  anthropomorphic  head  phantom  are  shown  in  hgure  7.  Figure  7(a)  shows  one 
slice  of  the  reconstructed  images  from  projection  data  acquired  with  the  10  mA  protocol. 
Figure  7(c)  shows  the  reconstructed  image  from  the  PWLS-processed  10  mA  sinogram. 
Figure  7(d)  shows  the  same  slice  of  the  image  reconstructed  from  the  sinogram  obtained  with 
80  mA.  It  can  be  observed  that  noise  in  10  mA  CT  images  is  efficiently  suppressed  after  the 
sinogram  is  processed  by  the  PWLS  algorithm.  The  processed  low-dose  CT  (10  mA)  image 
is  very  similar  to  its  corresponding  high-dose  image  (80  mA)  by  visual  judgment.  Standard 
deviation  of  the  noise  in  a  uniform  ROI  (as  indicated  by  an  arrow  in  hgure  7(d))  is  2.8  x  10“^in 
a  low-dose  (10  mA)  image  and  0.951  x  10“^  in  its  coiTesponding  high-dose  image.  After 
the  low-dose  sinogram  is  processed  by  the  PWLS  algorithm  with  a  smoothing  parameter 
P  —  0.05,  the  standard  deviation  of  the  same  ROI  is  0.955  x  10“^,  which  is  fairly  close  to  the 
noise  level  of  the  80  mA  image. 

To  further  illustrate  how  the  edge  information  is  affected  by  the  PWLS  sinogram 
smoothing,  in  hgure  7(e)  we  show  the  difference  image  between  hgures  7(a)  and  (c).  In 
the  difference  image,  random  noise  is  dominant  and  no  edge  or  structure  can  be  observed. 
This  indicates  that  the  edge  information  is  well  preserved  in  the  PWLS-processed  images. 

4.  Discussion 

Generally,  noise  reduction  for  CT  imaging  can  be  performed  in  three  spaces:  projection  data 
(either  before  or  after  logarithmic  transform),  hltered  projection  data  (before  backprojection 
operation  during  reconstruction)  and  reconstructed  CT  images.  During  hltering  and 
backprojection  operation,  the  noise  properties  will  change  signihcantly.  Then  noise  modeling, 
such  as  distribution  of  noise  and  variance  of  noise,  is  difficult  in  hltered  projection  data  and 
reconstructed  image.  Therefore,  in  this  work  we  chose  to  work  on  the  log-transformed  data 
to  fully  utilize  the  noise  model  of  the  projection  data  in  the  Radon  space  (Li  et  al  2004,  Wang 
et  al  2008). 
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Figure  6.  One  slice  of  the  FDK  reconstructed  image  of  the  CatPhan®  600  phantom  containing 
several  circular  objects  with  different  intensities:  (a)  from  projection  images  acquired  with  10  mA 
tube  current,  (b)  after  the  sinogram  acquired  with  10  mA  tube  current  is  processed  by  the  PWLS 
algorithm  and  (c)  from  projection  images  acquired  with  80  mA  tube  current. 


Accurate  noise  modeling  of  measurement  is  fundamentally  important  in  statistics-based 
image  processing  algorithms.  Meanwhile,  the  regularization  term  also  plays  an  important 
role  in  the  performance  of  the  algorithm.  In  CT  sinogram  processing,  a  commonly  used 
regularization  takes  a  quadratic  form  with  equal  weights  for  neighbors  of  an  equal  distance 
(La  Riviere  2005,  La  Riviere  and  Billmire  2005,  Li  et  al  2004,  Wang  et  al  2006).  Such  a 
quadratic  penalty  simply  encourages  the  equivalence  between  neighbors  without  considering 
discontinuities  in  the  image  and  may  lead  to  over-smoothing  around  sharp  edges  or  boundaries. 
In  the  presented  algorithm,  we  proposed  an  anisotropic  penalty  to  consider  the  difference 
among  neighbors.  The  idea  was  inspired  by  the  well-known  anisotropic  diffusion  filter  (Perona 
and  Malik  1990),  in  which  the  gradient  controls  the  strength  of  diffusion  among  neighbors.  The 
coupling  between  neighbors  should  be  smaller  if  the  absolute  value  of  difference  between  them 
is  smaller  and  this  kind  of  neighbors  should  contribute  less  to  the  solution  of  the  concerned 
pixel  (see  equation  (6)).  There  are  many  choices  that  satisfy  this  behavior  of  weighting.  In  this 
work,  the  form  of  the  anisotropic  weight  was  chosen  the  same  as  the  conduction  coefficients 
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Figure  7.  One  slice  of  FDK  reconstructed  image  of  the  anthropomorphic  head  phantom:  (a) 
from  projection  images  acquired  with  10  mA  tube  cuiTent,  (b)  using  a  low-pass  Hanning  filter 
with  cutoff  80%  Nyquist  frequency,  (c)  after  the  sinogram  acquired  with  10  mA  tube  current  is 
processed  by  the  PWLS  algorithm,  (d)  from  projection  images  acquired  with  80  mA  tube  current 
and  (e)  difference  image  between  (d)  and  (c). 


in  the  anisotropic  diffusion  filter  (Perona  and  Malik  1990).  By  such  a  choice,  the  anisotropic 
quadratic  form  penalty  discourages  the  equivalence  between  neighbors  if  the  gradient  between 
them  is  large,  and  the  edges  or  boundaries  in  the  image  will  be  better  preserved.  This  effect  is 
similar  to  that  of  anisotropic  coefficients  in  the  diffusion  filter. 

In  the  presented  method,  the  reconstruction  of  CT  images  was  performed  by  an  analytical 
FDK  algorithm  for  its  speed  and  accuracy.  During  the  FDK  reconstruction  process,  noise 
can  also  be  suppressed  by  using  a  low-pass  hlter.  It  has  been  reported  (Li  et  al  2004,  La 
Riviere  2005)  that  a  statistics-based  sinogram  smoothing  algorithm  plus  FBP  reconstruction 
is  superior  to  conventional  low-pass  filters  for  noise  suppression  of  2D  fan-beam  CT.  In  this 
work,  we  also  reconstructed  the  CT  image  of  the  anthropomorphic  head  phantom  using  a 
Hanning  filter  with  a  cutoff  at  80%  Nyquist  frequency,  see  figure  7(b).  It  can  be  observed  that 
the  image  reconstructed  from  the  PWLS-processed  sinogram  is  superior  to  the  result  of  the 
Hanning  hlter  in  terms  of  noise  suppression  and  structure  preservation. 

Similar  to  the  cutoff  frequency  in  the  conventional  low-pass  hlter  during  reconstruction, 
there  is  also  a  free  parameter  p  in  the  presented  method  which  controls  the  trade-off  of  the 
noise  level  and  the  structure  preservation  in  reconstructed  images.  In  this  work,  the  choice 
of  p  is  determined  by  the  visual  judgment.  The  optimal  choice  of  the  parameter  P  can  be 
determined  by  more  sophisticated  ways  such  as  the  received  operating  characteristic  (ROC) 
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Study.  Nevertheless,  the  parameter  /0  can  be  chosen  according  to  the  noise  level  of  the  sinogram 
because  from  equation  (6), the  solution  for  the  ideal  sinogram,  we  can  see  that  the  parameter 
P  and  variance  erf  are  always  coupled  together.  The  noise  level  of  the  projection  data  is 
mainly  determined  by  two  factors:  incident  photon  number  and  line  integrals.  As  such,  the 
parameter  p  could  be  optimized  at  a  certain  mAs  level  and  treatment  site  for  patients  of  a 
similar  size.  In  this  work,  however,  the  parameter  is  chosen  empirically,  which  is  justihable 
when  the  dependence  of  the  parameter  on  the  noise  level  is  weak. 

The  method  presented  in  this  paper  is  based  on  the  noise  properties  of  the  sinogram, 
and  the  smoothing  constraint  or  penalty  is  applied  to  the  sinogram  domain.  Based  on  the 
same  noise  model,  the  smoothing  constraint  can  also  be  applied  to  the  CT  image  domain,  and 
the  statistical  iterative  reconstruction  (SIR)  algorithm  can  be  used  to  obtain  the  attenuation 
coefficient  map  by  minimizing  the  objective  function.  The  SIR-based  algorithms  showed 
some  advantages  over  the  conventional  filtered  backprojection  method  for  multi-slice  helical 
CT  (Thibault  et  al  2007).  However,  an  obstacle  for  practical  use  of  SIR  is  the  computation 
burden,  especially  for  large  volume  datasets  of  CBCT.  It  takes  more  than  10  h  to  reconstruct 
the  typical  volume  of  multi-slice  helical  CT  using  SIR  (Thibault  et  al  2007).  It  takes  only 
3  s  for  the  presented  sinogram  smoothing  method  to  process  one  projection  image  using 
a  PC  with  3.0  GHz  CPU.  Parallel  computing  can  speed  up  both  SIR  and  sinogram-based 
algorithms  significantly  using  the  cell  broadband  engine  (Knaup  et  al  2006)  and  PC  cluster 
(Li  et  al  2005b).  It  is  possible  to  achieve  clinically  acceptable  time  for  the  presented  sinogram 
smoothing  algorithm  through  parallel  computation.  It  is  an  interesting  research  topic  to 
quantitatively  compare  the  performance  of  the  SIR-based  CBCT  reconstruction  algorithm  and 
the  statistics-based  sinogram  smoothing  method. 

When  CBCT  is  used  for  patient  setup  and  target  localization  during  radiotherapy,  some 
extra  information  may  be  taken  into  account  for  dose  and  noise  reduction.  For  example,  a 
complete  CT  volume  dataset  with  high  clarity  used  for  treatment  planning  is  usually  available 
before  the  treatment.  This  will  provide  strong  a  priori  information  of  the  patient  before  each 
CBCT  scan.  Prior  information  of  planning  3D  CT  has  been  proved  useful  to  improve  the 
image  quality  of  4D  CBCT  (Li  et  al  2007).  It  is  expected  that  the  radiation  dose  of  CBCT 
used  for  radiotherapy  can  be  further  reduced  by  incorporating  the  planning  CT  information 
into  the  image  restoration  or  reconstruction  algorithms. 

In  the  report  of  the  AAPM  task  group  75  (Murphy  et  al  2007),  several  dose  reduction 
strategies  for  image-guided  radiotherapy  were  discussed.  For  CBCT,  dose  reduction  can  be 
achieved  by  narrowing  field  of  view  to  avoid  delivering  radiation  to  an  unnecessary  region  of 
the  patient  (Murphy  et  al  2007).  Compared  with  these  hardware-based  approaches,  software 
approaches  (such  as  the  one  proposed  in  this  paper)  provide  a  more  cost-effective  means 
for  dose  reduction  of  CBCT.  In  addition  to  the  statistics-based  reconstruction  and  restoration 
algorithms,  advanced  analytical  CBCT  reconstruction  algorithms  (Leng  et  al  2007,  Zhuang 
et  al  2006,  Zou  and  Pan  2004,  Zou  et  al  2005)  may  further  improve  the  low-dose  CBCT  image 
quality. 

5.  Conclusion 

A  PWLS  algorithm  with  non-uniform  weights  was  proposed  to  reduce  noise  in  low-dose 
onboard  CBCT.  In  this  method,  the  sinogram  was  first  processed  according  to  the  PWLS 
criterion.  The  weight  for  each  measurement  was  chosen  as  sinogram  datum  variance,  where 
variance  can  be  estimated  accurately  according  to  the  sinogram  noise  model.  To  preserve 
edge  information  during  noise  reduction,  we  proposed  an  anisotropic  quadratic  form  penalty. 
The  quadratic  form  penalty  encourages  equivalence  between  neighbors  and  the  anisotropic 
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penalty  provides  the  mechanism  to  control  the  influence  of  different  neighbors  according  to 
its  corresponding  gradient.  The  effectiveness  of  the  proposed  method  is  demonstrated  by 
two  experimental  phantom  studies.  The  quality  of  the  10  mA  CT  image  after  its  sinogram 
processed  by  the  PWLS  algorithm  is  comparable  to  the  image  obtained  with  the  80  mA 
protocol.  These  experimental  results  indicate  that  it  is  possible  to  reduce  the  CBCT  radiation 
dose  by  a  factor  of  1  /8  without  loss  of  useful  information  for  radiotherapy. 
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On-board  cone-beam  computed  tomography  (CBCT)  is  a  new  imaging  technique  for  radiation 
therapy  guidance,  which  provides  volumetric  information  of  a  patient  at  treatment  position.  CBCT 
improves  the  setup  accuracy  and  may  be  used  for  dose  reconstruction.  However,  there  is  great 
concern  that  the  repeated  use  of  CBCT  during  a  treatment  course  delivers  too  much  of  an  extra  dose 
to  the  patient.  To  reduce  the  CBCT  dose,  one  needs  to  lower  the  total  mAs  of  the  x-ray  tube  current, 
which  usually  leads  to  reduced  image  quality.  Our  goal  of  this  work  is  to  develop  an  effective 
method  that  enables  one  to  achieve  a  clinically  acceptable  CBCT  image  with  as  low  as  possible 
mAs  without  compromising  quality.  An  iterative  image  reconstruction  algorithm  based  on  a  penal¬ 
ized  weighted  least-squares  (PWLS)  principle  was  developed  for  this  purpose.  To  preserve  edges  in 
the  reconstructed  images,  we  designed  an  anisotropic  penalty  term  of  a  quadratic  form.  The  algo¬ 
rithm  was  evaluated  with  a  CT  quality  assurance  phantom  and  an  anthropomorphic  head  phantom. 
Compared  with  conventional  isotropic  penalty,  the  PWLS  image  reconstruction  algorithm  with 
anisotropic  penalty  shows  better  resolution  preservation.  ©  2009  American  Association  of  Physi¬ 
cists  in  Medicine.  [DOI;  10.1118/1.3036112] 

Key  words;  cone-beam  CT,  low-dose,  iterative  reconstruction,  PWLS,  edge-preserving  penalty 


I.  INTRODUCTION 

Integration  of  the  cone-beam  computed  tomography  (CBCT) 
with  a  linear  accelerator'  makes  it  possible  to  acquire  volu¬ 
metric  image  of  high  spatial  resolution  for  a  patient  at  treat¬ 
ment  position.  There  is  growing  interest  in  using  on-board 
CBCT  for  patient  setup  and  dose  reconstruction.  Repeated 
use  of  CBCT  during  a  treatment  course  has  raised  concern  of 
the  extra  radiation  dose  delivered  to  patients.^’"'  One  cost- 
effective  way  to  reduce  the  CBCT  dose  is  to  acquire  CT  with 
a  lower  mAs  protocol.  However,  image  quality  will  degrade 
dramatically  due  to  excessive  noise,^’^  rendering  the  low- 
mAs  CBCT  a  less  attractive  option  for  the  therapeutic  guid¬ 
ance. 

In  this  work,  we  incorporate  the  noise  properties  of  CBCT 
log-transformed  projection  data'^’^  in  a  statistical  iterative  im¬ 
age  reconstruction  algorithm  to  improve  the  low-dose  CBCT 
image  quality.  Compared  with  analytical  reconstruction  algo¬ 
rithms,  a  major  advantage  of  iterative  algorithms  is  that  it 
takes  into  consideration  the  image  physics,  noise  properties, 
and  imaging  geometry  elegantly.  Advantages  of  iterative  re¬ 
construction  algorithms  have  been  demonstrated  in  the  image 
reconstruction  of  emission  tomographic  images.^ How¬ 
ever,  when  applying  iterative  reconstruction  algorithms  for 
CT  imaging,  '  long  computational  time  may  pose  a  chal¬ 
lenge  for  their  clinical  applications.  With  the  development  of 
fast  computers  and  dedicated  hardwares,'*’*^  iterative  recon¬ 
struction  algorithms  may  be  used  for  clinical  CT  reconstruc¬ 


tion  in  the  near  future.  Recently,  iterative  image  reconstruc¬ 
tion  algorithms  have  demonstrated  superior  performance  for 
reconstruction  of  the  multislice  helical  CT  (Ref.  20)  and  car- 
diac  micro-CT.  Prototype  products  based  on  iterative  re¬ 
construction  methods  have  been  exhibited  by  major  CT  ven¬ 
dors  in  a  number  of  national  and  international  meetings. 

In  statistics-based  iterative  reconstruction  algorithms,  to¬ 
mographic  images  are  reconstructed  by  minimizing  or  maxi¬ 
mizing  a  cost  function,  which  is  constructed  based  on  noise 
characteristics  of  the  measured  data.  There  are  usually  two 
terms  in  the  objective  function.  The  first  term  models  the 

statistics  of  measured  data  and  the  second  term  reflects  a 

1  22 

prior  information  to  regularize  the  solution.  Many  efforts  ’ 
have  been  devoted  to  investigate  the  noise  models  of  the 
measurements  in  CT.  Accurate  noise  modeling  is  a  prerequi¬ 
site  of  a  statistical  iterative  reconstruction  algorithm.  The 
second  term,  i.e.,  the  regularization  term,  also  plays  an  im¬ 
portant  role  for  successful  image  reconstruction.  One  com¬ 
mon  choice  of  the  regularization  term  is  the  Gaussian  Mar¬ 
kov  random  field  in  quadratic  form.  Such  quadratic 

penalty  with  equal  weights  for  neighbors  of  equal  distance 
encourages  equivalence  between  neighbors  without  consid¬ 
ering  discontinuities  in  images,  which  may  lead  to  over¬ 
smoothing  around  edges  or  boundaries.  Several  edge¬ 
preserving  regularization  methods  have  been  proposed  to 
address  this  problem.  For  example,  the  edge-preserving  Hu¬ 
ber  penalty,  which  penalizes  neighbors  of  small  differences 
quadratically  while  applying  a  linear  penalty  on  neighbors  of 
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larger  difference,  has  been  used  by  Elbakri  and  Fessler^^  for 
CT  image  reconstruction  and  by  Chlewicki  et  al.  for  posi¬ 
tron  emission  tomography  reconstruction.  A  line  process  has 
been  introduced  by  Geman  and  Geman  to  define  the  edge 
lattice  during  Bayesian  restoration  of  images.  Geman  and 
Reynolds  proposed  a  finite  asymptotic  edge-preserving 
function  and  Charbonnier  et  introduced  an  auxiliary 
variable  in  the  prior  constraint  to  mark  the  discontinuities  in 
the  images.  These  modifications  make  the  objective  function 
nonquadratic  and  complicate  the  computation.  In  this  work, 
we  propose  a  quadratic  regularization  term  with  anisotropic 
weights  for  different  neighbors.  The  role  of  the  anisotropic 
penalty  is  to  discourage  the  equivalence  between  neighbors 
if  the  gradient  is  large;  thus  the  edges  or  discontinuities  will 
be  better  preserved  in  the  final  reconstructed  image. 

In  the  following  sections,  we  first  introduce  the  penalized 
weighted  least  squares  (PWLS)  objective  function  for  image 
reconstruction  of  CBCT  based  on  the  noise  properties  of 
CBCT  projection  data.  We  then  describe  the  proposed  aniso¬ 
tropic  penalty  in  details.  In  Sec.  Ill,  the  evaluation  of  the 
proposed  algorithm  is  presented  using  a  quality  assurance 
phantom  and  an  anthropomorphic  head  phantom,  followed 
by  the  discussion  in  Sec.  IV  and  the  conclusion  in  Sec.  V. 

II.  METHODS  AND  MATERIALS 
II.A.  PWLS  image  reconstruction 

Noise  in  x-ray  CT  projection  data  after  logarithm  trans¬ 
form  follows  approximately  Gaussian  distribution  and  the 
variance  of  the  noise  can  be  determined  by  an  exponential 
formula^’* 

=  exp(p,.)/V,o,  (1) 

where  V,o  is  the  incident  photon  number  at  detector  bin  i,  p, 
and  (jf  is  the  mean  and  variance  of  projection  datum 
respectively.  Based  on  the  noise  properties  of  CT  projection 
data,  the  PWLS  cost  function  in  the  image  domain  can  be 
written  as 

^{lj)  =  {p-Ap)'l.~^{p-Ap)  +  l3R{p).  (2) 

The  first  term  in  Eq.  (2)  is  the  weighted  least-squares  crite¬ 
rion,  where  p  is  the  vector  of  log-transformed  projection 
measurements,  and  pu  is  the  vector  of  attenuation  coefficients 
to  be  reconstructed.  Operator  A  represents  the  system  or  pro¬ 
jection  matrix.  The  element  of  is  the  length  of  the  inter¬ 
section  of  projection  ray  i  with  pixel  j  and  it  is  calculated  by 

'll 

a  fast  ray-tracing  technique.  In  our  implementation,  the  el¬ 
ement  of  matrix  A  was  precomputed,  stored  as  a  file,  and 
used  as  a  lookup  table  later.  The  projection  data  p  and  the 
attenuation  map  p,  is  related  by  p=Api.  X  is  a  diagonal  ma¬ 
trix  with  the  ;th  element  of  of,  i.e.,  an  estimate  of  the  vari¬ 
ance  of  measured  y,  at  detector  bin  i  which  can  be  estimated 
from  the  measured  projection  data  according  to  Eq.  (1).  The 
element  of  the  diagonal  matrix  plays  the  role  of  weighting  in 
the  WLS  cost  function  and  it  determines  contribution  of  each 
measurement.  The  symbol  '  denotes  the  transpose  operator. 
The  second  term  in  Eq.  (2)  is  a  smoothness  penalty  or  a  prior 
constraint,  where  (3  is  the  smoothing  or  penalty  parameter 


which  controls  relative  contribution  from  the  measurement 
and  prior  constraint.  The  image  reconstruction  task  is  to  find 
attenuation  map  p,  by  minimizing  the  objective  function  (2) 
with  a  positive  constraint.  The  Gaussian-Seidel  updating 
strategy  was  used  for  the  minimization  and  details  about  the 
implementation  are  described  in  the  Appendix. 


II. B.  Edge-preserving  anisotropic  penaity 

The  prior  constraint  in  Eq.  (2)  enforces  a  roughness  pen¬ 
alty  on  the  solution.  The  quadratic  penalty  with  equal 
weights  for  neighbors  of  the  same  distance  has  been  used 

16  17  23  25 

widely  for  iterative  image  reconstruction  ’  ’  ’ 

R{p)  =  p'Rp.=  -^  2  Wj„JyPj-pS,  (3) 

^  j  m&Nj 


where  index  j  runs  over  all  image  elements  in  the  image 
domain,  Nj  represents  the  set  of  neighbors  of  the  y'th  image 
pixel.  The  parameter  Wj,„  was  set  to  1  for  first-order  neigh¬ 
bors  and  1/^2  for  second-order  neighbors  in  previous 
applications.  ’  ’  This  type  of  penalty  only  takes  distance 
information  of  the  neighbors  into  account.  That  is,  the  neigh¬ 
bors  of  the  same  distance  play  an  equivalent  role  in  regular¬ 
izing  the  solution,  and  vise  versa.  A  major  shortcoming  of 
the  approach  is  that  the  regularization  does  not  take  the  dif¬ 
ference  in  intensities  of  the  neighboring  voxles  (e.g.,  edges 
or  discontinuities)  into  account,  which  may  lead  to  an  over¬ 
smoothed  solution  for  reconstructed  images.  To  overcome 
this  limitation,  we  propose  an  anisotropic  penalty  to  regular¬ 
ize  the  solution.  In  this  formulation,  the  weight  is  smaller  if 
the  difference  between  a  neighbor  and  the  concerned  voxel  is 
larger,  since  the  coupling  between  two  such  neighbors  is 
smaller.  There  are  many  choices  that  satisfy  this  behavior  of 
weighting.  In  this  work,  we  chose  the  form  of  Wj„,  to  be  the 
same  as  the  conduction  coefficient  in  the  well-known  aniso- 


32 

tropic  diffusion  filter.  The  weight  Wj„j  can  be  written  as 


-  P-n 


(4) 


where  the  gradient  and  the  parameter  S  determine  the 
strength  of  the  diffusion  during  each  iteration.  The  parameter 
S  can  be  set  either  by  hand  or  to  the  value  at  90%  of  the 
histogram  of  the  gradient  magnitude  of  the  image  to  be  pro¬ 
cessed.  In  this  work,  we  set  the  value  of  S  to  be  90%  of  the 
histogram  of  the  gradient  magnitude  of  the  FDK  recon¬ 
structed  image  (which  is  used  as  the  initial  during  iterative 
reconstruction). 


II.C.  CBCT  data  acquisition 

Cone-beam  CT  projection  data  were  acquired  by  an  Acu¬ 
ity  simulator  (Varian  Medical  Systems,  Palo  Alto,  CA).  The 
number  of  projections  for  a  full  360°  rotation  is  680  and  the 
total  time  for  the  acquisition  of  one  full  circle  of  the  projec¬ 
tion  data  is  about  1  min.  The  dimension  of  each  acquired 
projection  image  is  397  mm  X  298  mm,  containing  1024 
X  768  pixels.  To  save  computational  time  during  iterative  re¬ 
construction,  the  projection  data  at  each  projection  view 
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Fig.  1.  The  bead  point  object  in  the  CatPhan®  600  phantom  was  used  to 
calculate  the  MTF  of  reconstructed  images.  Display  window: 
[0,0.03]  mm“'.  The  white  square  in  the  image  indicates  the  region  used  to 
calculate  the  standard  deviation. 


were  downsampled  by  a  factor  of  2  and  only  16  out  of  768 
projection  data  along  the  axial  direction  were  chosen  for  re¬ 
construction.  The  system  has  a  field  of  view  of  25  cm 
X  25  cm  (full-fan  mode)  in  the  transverse  plane  and  17  cm  in 
the  longitudinal  direction,  which  can  be  increased  to  45  cm 
X  45  cm  in  the  transverse  plane  by  shifting  the  detector  lat¬ 
erally  (half-fan  mode). 

Two  phantoms  were  used  to  evaluate  the  performance  of 
the  proposed  PWLS  algorithm.  The  first  is  a  commercial 
calibration  phantom  CatPhan®  600  (The  Phantom  Labora¬ 
tory,  Inc.,  Salem,  NY).  The  second  is  an  anthropomorphic 
head  phantom.  In  both  phantom  studies,  the  tube  voltage  was 
set  to  125  kVp.  The  x-ray  tube  current  was  set  at  10  mA  and 
the  duration  of  the  x-ray  pulse  at  each  projection  view  was 
10  ms  during  the  acquisition  of  low-dose  CBCT  projection 
data.  During  acquisition  of  the  corresponding  high-dose 
CBCT  image,  the  tube  current  was  set  at  80  mA  and  the 
duration  of  the  x-ray  pulse  was  set  at  12  ms.  The  projection 
data  were  acquired  in  full-fan  mode  and  the  full-fan  bow-tie 
filter  was  used  for  both  phantoms.  The  distance  of  source-to- 
axis  is  100  cm  and  source-to-detector  distance  of  150  cm. 
The  size  of  reconstructed  image  is  350  X  350  X  16  and  voxel 
size  is  0.776X0.776X0.776  mm^. 

II. D.  Performance  evaluation 

We  used  the  CatPhan®  600  phantom  to  study  the  spatial 
resolution  of  images  reconstructed  by  different  algorithms. 
The  CTP591  module  of  the  CatPhan®  600  phantom  contains 
a  bead  point  object  with  a  diameter  of  0.28  mm  (see  Fig.  1). 
The  point  object  can  be  used  to  calculate  the  modulation 
transfer  function  (MTF)  which  characterizes  the  spatial  res¬ 
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Fig.  2.  MTF  curves  of  two  PWLS  iterative  image  reconstruction  algorithms 
with  different  smoothing  parameters.  Curves  in  the  blue  color  are  results  of 
reconstruction  using  an  isotropic  quadratic  penalty.  Curves  in  the  black  color 
are  results  of  reconstruction  using  the  edge-preserving  anisotropic  quadratic 
penalty. 

olution  of  images.  The  reconstructed  image  contains  the 
point  object  provides  the  point-spread-function  for  each  re¬ 
construction  algorithm  and  the  in-plane  MTF  can  be  ob¬ 
tained  by  calculating  two-dimension  Fourier  transform  and 
averaging  over  27r  angles.  A  lOX  10  matrix  centered  about 
the  point  object  was  used  to  calculate  the  MTF  after  the 
background  value  (which  can  be  estimated  by  averaging  the 
values  of  a  uniform  background  region)  was  subtracted  from 
the  value  of  each  pixel. 

In  the  CTP404  module  of  the  CatPhan®  600,  there  are 
several  circles  of  different  intensities  which  can  be  used  to 
quantify  the  contrast-to-noise  (CNR)  of  the  reconstructed  im¬ 
ages  in  different  reconstructions.  We  selected  a  low-contrast 
region  of  interest  (ROI)  for  calculation  of  the  CNR  in  the 
image  reconstructed  by  different  algorithms  since  a  low- 
contrast  region  is  of  interest  in  CT  imaging.  The  contrast  was 
calculated  as  the  absolute  difference  between  the  mean  value 
of  the  region  inside  the  ROI  and  the  mean  value  of  the  uni¬ 
form  background  region.  The  noise  level  was  characterized 
by  the  standard  deviation  of  a  uniform  area  of  size  15  pixels 
by  15  pixels.  The  CNR  was  defined  as  the  contrast  divided 
by  the  standard  deviation. 

III.  RESULTS 

III.A.  CatPhan®  600  phantom 
III. A.  1.  MTF  measurement 

Figure  2  shows  the  MTFs  of  two  iterative  reconstruction 
algorithms  with  different  smoothing  parameters  (3  ranging 
from  1.0  X  10"^  to  30  X  10"^.  It  can  be  observed  that  the  spatial 
resolution  of  the  reconstructed  image  using  an  isotropic  qua¬ 
dratic  penalty  decreases  as  smoothing  strength  increases.  The 
frequency  of  50%  MTF  for  the  iterative  reconstruction  using 
the  isotropic  penalty  decreases  from  5.25  to  2.91  1/cm  as 
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(e)  (f) 


Fig.  3.  CBCT  of  the  CatPhan®  600  phantom:  (a)  analytical  FDK  recon¬ 
structed  image  from  projection  data  acquired  using  a  low-dose  protocol 
(10  mA/10  ms)  and  (b)  a  high-dose  protocol  (80  mA/12  ms);  (c)  PWLS 
iterative  image  reconstruction  with  the  isotropic  quadratic  penalty  from  pro¬ 
jection  data  acquired  using  a  low-dose  protocol  and  (d)  with  the  proposed 
anisotropic  penalty;  (e)  PWLS  iterative  image  reconstruction  with  the  Huber 
penalty;  (f)  analytical  FDK  reconstructed  image  after  low-dose  projection 
data  processed  by  the  PWLS  criterion  (Ref.  34)  with  a  smoothing  parameter 
of  0.09.  Display  window:  [0,0.02]  mm“^ 

the  smoothing  parameter  increases  from  1.0X10‘*to30 
X  10^.  In  contrast,  MTF  curves  of  the  image  reconstructed 
using  the  proposed  edge-preserving  anisotropic  quadratic 
penalty  are  clustered  together  with  various  smoothing  param¬ 
eters.  This  indicates  that  the  spatial  resolution  in  images  re¬ 
constructed  using  the  anisotropic  quadratic  penalty  is  better- 
preserved. 

III.A.2.  Full  width  at  half  maximum  measurement 

We  then  tested  the  proposed  algorithm  on  the  CTP404 
module  of  the  CatPhan®  600  phantom.  A  representative  slice 
of  the  CBCT  images  obtained  by  different  reconstruction 
methods  are  shown  in  Fig.  3.  Figure  3(a)  is  the  low-dose 
image  reconstructed  by  analytical  FDK  algorithm.  It  can  be 
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observed  that  noise  level  is  high  in  this  low-dose  CBCT 
image.  Figure  3(c)  shows  the  image  reconstructed  by  the 
PWLS  algorithm  using  the  isotropic  quadratic  penalty  with 
the  penalty  parameter  (3=30  X  10^.  Figure  3(d)  displays  the 
image  reconstructed  by  the  PWLS  algorithm  using  the  pro¬ 
posed  edge-preserving  anisotropic  penalty  with  the  same 
penalty  parameter.  The  noise  in  the  images  reconstructed  by 
iterative  algorithms  is  greatly  suppressed  compared  with  the 
image  reconstructed  using  the  analytical  method.  It  is  seen 
that  the  edges  were  blurred  in  Fig.  3(c),  as  indicated  by  the 
arrows  in  the  image.  This  is  not  surprising  since  an  isotropic 
quadratic  form  penalty  simply  encourages  equivalence 
among  neighbors  along  all  directions  without  considering  the 
boundary  information  presented  in  the  image.  However, 
edges  were  well  preserved  when  the  anisotropic  penalty  was 
used  as  the  constraint  (see  the  corresponding  area  in  Fig.  3(d) 
indicated  by  the  arrows). 

To  quantitatively  analyze  the  gain  by  using  the  anisotropic 
penalty  in  the  iterative  reconstruction  algorithm,  we  then 
studied  the  full-width-at-half-maximum  (FWHM)  of  two 
pointlike  objects  (one  is  brighter  than  background  and  the 
other  one  is  darker  than  the  background)  in  the  reconstructed 
images.  Figure  4  shows  the  prohles  passing  through  two 
pointlike  objects  in  Figs.  3(c)  and  3(d).  Through  those  pro¬ 
hles,  it  can  be  observed  that  the  major  difference  between  the 
solutions  using  isotropic  and  anisotropic  penalties  is  nearby 
edges.  The  intensity  values  in  both  images  at  a  uniform  re¬ 
gion  are  nearly  identical;  see  line  along  value  0.016  in  both 
Figs.  4(a)  and  4(b).  It  can  also  be  observed  in  Fig.  4  that  the 
peak  of  the  prohle  from  the  image  reconstructed  using  the 
isotropic  quadratic  penalty  is  lower  than  that  from  the  recon¬ 
structed  image  using  the  anisotropic  penalty,  while  the  bot¬ 
tom  of  the  prohle  from  the  image  reconstructed  using  the 
isotropic  quadratic  penalty  is  higher  than  that  from  the  re¬ 
constructed  image  using  anisotropic  penalty.  These  observa¬ 
tions  show  that  there  is  a  signal  loss  when  an  image  is  re¬ 
constructed  by  the  PWLS  algorithm  using  the  isotropic 
penalty.  The  standard  deviation  of  a  uniform  region  (indi¬ 
cated  by  a  white  square)  is  0.50X10'^  in  Fig.  3(c)  and 
0.54  X  10'^  in  Fig.  3(d).  We  then  htted  the  prohle  to  a  Gauss¬ 
ian  functional.  The  FWHM  of  brighter  source  is  3.48  pixels 
for  the  image  reconstructed  with  the  isotropic  penalty  and 
3.19  pixels  for  the  anisotropic  penalty.  The  FWHM  for  the 
darker  source  is  3.71  pixels  for  the  image  reconstructed  with 
the  isotropic  penalty  and  3.48  pixels  for  the  anisotropic  pen¬ 
alty.  In  both  cases,  better  edge  preserving  was  observed  in 
the  image  reconstructed  using  the  anisotropic  penalty. 

III.A.3.  CNR  measurement 

Table  I  lists  the  CNR  of  two  iterative  reconstruction  algo¬ 
rithms  with  different  smoothing  parameters.  It  can  be  ob¬ 
served  that  in  both  reconstruction  algorithms  the  CNR  in¬ 
crease  as  smoothing  parameter  increases.  The  CNR  of  the 
image  reconstructed  using  the  anisotropic  penalty  is  slightly 
larger  than  that  of  the  isotropic  penalty  when  the  same 
smoothing  parameters  are  used.  However,  at  the  matched 
resolution  between  the  two  methods,  CNR  was  increased 
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O  Isotropic  FWHM  =  3.48  A  PDK  FWHM  =  3.29 
*  Anisotropic  FWHM  =3.19 


(a) 

o  Isotropic  FWHM  =  3.71  A  fDK  FWHM  =  3.50 
*  Anisotropic  FWHM  =  3.48 


(b) 
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Table  I.  CNRs  of  the  low-contrast  ROI  in  Fig.  3. 


yS  (XlOb 

1 

2.5 

7.5 

10 

20 

30 

PWLS  isotropic 

0.84 

0.98 

1.40 

1.60 

2.35 

3.01 

PWLS  anisotropic 

0.77 

0.82 

1.04 

1.20 

2.22 

2.83 

FDK  (10  mA/10  ms) 

0.95 

FDK  (80  mA/12  ms) 

2.66 

H{t)  = 


f-n, 

6»(|f|  -e)  +  ^/2, 


\t\<e 


(5) 


The  Huber  function  penalizes  the  difference  between  neigh¬ 
boring  pixels  quadratically  if  the  absolute  difference  pixel 
value  |f|  is  smaller  than  some  threshold  6  and  it  will  apply  a 
linear  penalty  to  the  larger  differences  of  \t\^  6  which  usu¬ 
ally  occur  at  edges. 

Figure  3(e)  shows  the  PWLS  reconstructed  CatPhan®  600 
phantom  by  using  the  Huber  penalty  with  threshold  6 
=  0.001  and  the  penalty  parameter  yS=35  X  10"^.  It  can  be  ob¬ 
served  that  the  edges  are  better  preserved  in  the  images  re¬ 
constructed  using  the  Huber  penalty  than  the  images  recon¬ 
structed  by  using  the  isotropic  quadratic  penalty.  To 
quantitatively  compare  the  performance  of  the  Huber  penalty 
and  the  anisotropic  quadratic  penalty,  we  calculated  the  MTF 
of  the  CTP591  module  of  the  CatPhan®  600  phantom  at 
matched  noise  level.  Figure  5  shows  the  MTF  curves  from 
the  proposed  anisotropic  quadratic  penalty  with  penalty  pa¬ 
rameter  yS=30X10'*  and  the  Huber  penalty  with  threshold 
0=0.001.  The  penalty  parameter  (3  in  the  PWLS  reconstruc¬ 
tion  with  Huber  penalty  was  set  at  35  X  10"^  so  that  the  noise 
level  in  the  reconstructed  image  is  matched  to  the  anisotropic 
quadratic  penalty.  The  standard  deviation  of  the  uniform  re¬ 
gion  (indicated  by  a  white  square  in  Fig.  1)  in  the  image 
reconstructed  using  the  Huber  penalty  is  0.73  X  10“^  and  is 
0.70  X  10“^  in  the  image  reconstructed  using  the  anisotropic 
quadratic  penalty.  MTF  curves  in  Fig.  5  show  that  the  aniso- 


Fig.  4.  Vertical  profile  through  column  139  in  Figs.  3(c)  and  3(d).  (a)  shows 
the  profile  through  the  hot  spot  and  (b)  shows  the  profile  through  the  cold 
spot.  Edges  are  better  preserved  by  using  the  anisotropic  penalty  as  mea¬ 
sured  by  the  FWHM. 


from  0.84  in  the  image  reconstructed  using  the  isotropic  pen¬ 
alty  with  ;S=1.0X10'*  to  2.83  in  the  image  reconstructed 
using  the  anisotropic  penalty  with  /8=30X  lO'*. 


III.A.4.  Comparison  study  with  the  Huber  penaity 

In  this  section,  we  compared  the  proposed  anisotropic 
quadratic  penalty  with  a  representative  edge-preserving  non¬ 
quadratic  penalty;  the  Huber  penalty.^®’^^  The  Huber  penalty 
function  has  the  following  form; 


Fig.  5.  MTF  curves  of  the  low-dose  CBCT  image  reconstmcted  by  different 
algorithms  at  a  match  noise  level. 
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Fig.  6.  MTF  curves  of  the  high-dose  CBCT  image  reconstructed  by  differ¬ 
ent  algorithms. 


tropic  quadratic  penalty  produces  better  image  resolution  at 
the  matched  noise  level.  The  advantage  of  the  anisotropic 
quadratic  penalty  may  be  attributed  to  that  the  Huber  func¬ 
tion  depends  on  a  hard  threshold  of  the  gradient  while  the 
anisotropic  quadratic  penalty  considers  the  gradient  informa¬ 
tion  continuously  by  introducing  Eq.  (4). 


III.A.5.  Comparison  with  high-dose  CBCT 

For  the  projection  data  acquired  with  a  tube  current  of 
80  mA  and  x-ray  pulse  duration  of  12  ms  protocol,  we  re¬ 
constructed  the  CBCT  image  using  the  analytical  FDK  algo¬ 
rithm.  We  hrst  compared  the  MTF  of  the  image  recon¬ 
structed  by  the  analytical  FDK  algorithm  with  the  image 
reconstructed  by  iterative  PWLS  algorithms  at  a  matched 
noise  level.  The  standard  deviation  of  the  uniform  area  is 
5.95  X  10“'*  in  the  FDK-reconstructed  image.  By  setting  the 
smoothing  parameter  ;S=2X10'*,  the  standard  deviation  of 
the  same  region  in  the  PWLS -reconstructed  image  is  5.96 
X  10“'*  with  the  isotropic  penalty  and  5.98  X  10“'*  with  the 
anisotropic  penalty.  Figure  6  shows  the  MTF  curves  from  the 
image  reconstructed  by  FDK  and  the  iterative  PWLS  algo¬ 
rithms.  It  can  be  observed  that  the  MTF  of  the  PWLS  algo¬ 
rithm  with  the  anisotropic  penalty  is  slightly  better  than  that 
of  the  FDK  algorithm,  whereas  the  MTF  of  the  FDK  algo¬ 
rithm  is  better  than  that  of  the  PWLS  algorithm  with  isotro¬ 
pic  penalty.  This  demonstrates  that  better  spatial  resolution  is 
achieved  by  the  PWLS  algorithm  using  the  proposed  edge¬ 
preserving  anisotropic  penalty.  The  same  trend  can  also  be 
seen  from  the  prohles  through  the  pointlike  objects  in  the 
CTP404  module  (Fig.  4).  The  FWHM  obtained  from  the  ht- 
ted  Gaussian  function  also  shows  that  better  spatial  reso¬ 
lution  is  achieved  by  using  the  PWLS  image  reconstruction 
algorithm  with  the  anisotropic  penalty. 

From  Table  I,  it  is  seen  that  the  CNR  of  the  PWLS  recon¬ 
structed  low-dose  image  using  the  anisotropic  penalty  with 
the  penalty  parameter  yS=30X  10'*  is  2.83,  which  is  slightly 
higher  than  that  of  the  FDK  reconstructed  high-dose 
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image — 2.66.  In  Fig.  6,  we  also  show  the  MTF  curve  of  the 
PWLS  algorithm  using  the  anisotropic  penalty  with  the  pen¬ 
alty  parameter  f3=30X  10'*.  It  can  be  observed  that  the  MTF 
of  the  PWLS  algorithm  with  the  anisotropic  penalty  is  com¬ 
parable  to  that  of  the  FDK  algorithm.  This  result  suggests 
that  the  PWLS  iterative  image  reconstruction  with  the  aniso¬ 
tropic  penalty  is  capable  of  producing  images  with  a  CNR 
comparable  to  FDK-reconstructed  high-dose  images  using 
only  about  1/10  dose  without  sacrihcing  image  spatial  reso¬ 
lution. 


III.B.  Anthropomorphic  head  phantom 

Results  of  the  anthropomorphic  head  phantom  are  shown 
in  Fig.  7.  Figure  7(a)  shows  one  slice  of  the  image  recon¬ 
structed  by  the  analytical  FDK  algorithm  from  projection 
data  acquired  using  a  low-dose  protocol  (10  mA/ 10  ms); 
Fig.  7(b)  is  the  FDK  reconstructed  image  for  the  same  phan¬ 
tom  acquired  with  a  high-dose  protocol  (80  mA/ 12  ms). 
Figure  7(c)  shows  the  same  slice  of  a  low-dose  CBCT  image 
reconstructed  by  the  PWLS  iterative  algorithm  using  the  iso¬ 
tropic  penalty  with  the  penalty  parameter  of  (3=30  X  10'*  and 
Fig.  7(d)  shows  the  low-dose  CBCT  image  reconstructed  by 
the  PWLS  image  reconstruction  algorithm  using  the  aniso¬ 
tropic  penalty  with  the  same  penalty  parameter.  Figure  7(e) 
shows  the  PWLS  reconstructed  low-dose  CBCT  image  using 
the  edge-preserving  Huber  penalty  with  threshold  0=0.001 
and  the  penalty  parameter  (3=35  X  10'*.  It  can  be  observed 
that  noise  in  low-dose  CT  images  is  efficiently  suppressed  in 
images  reconstructed  by  the  PWLS  iterative  reconstruction 
algorithms.  The  quality  of  low-dose  CBCT  reconstructed  by 
the  PWLS  with  anisotropic  penalty  is  comparable  to  that  of 
the  high-dose  FDK  reconstructed  image.  With  the  anisotropic 
penalty  in  the  PWLS  iterative  reconstruction,  edges  are  bet¬ 
ter  preserved  in  the  reconstructed  image.  In  the  regions  indi¬ 
cated  by  arrows  in  Fig.  7(d),  it  is  seen  that  the  structure  is 
well  preserved  in  the  image  reconstructed  using  the  aniso- 
tropically  penalized  PWLS  algorithm.  The  structure  is 
blurred  if  the  isotropic  penalty  was  used  during  the  PWLS 
reconstruction.  This  observation  is  consistent  with  the  quan¬ 
titative  evaluation  using  the  CatPhan®  600  phantom. 


IV.  DISCUSSION 

The  weighted  least-squares  criterion  reflects  that  the  mea¬ 
sured  data  with  a  lower  SNR  will  contribute  less  to  the  esti¬ 
mation  of  the  attenuation  map.  The  PWLS  objective  function 
is  equivalent  to  the  penalized  maximum  likelihood  or  maxi¬ 
mum  a  posteriori  criterion  for  Gaussian  distributed  noise. 
This  is  consistent  with  the  observations  from  repeated  mea¬ 
surements  of  x-ray  CT  projection  data  after  logarithm 
transform.  ’  The  PWLS  criterion  for  the  CT  projection  data 
can  also  be  derived  from  Poisson  noise  model  of  detector 
counts  using  the  second-order  Taylor  expansion.  *®  Measure¬ 
ment  of  x-ray  counts  can  be  modeled  more  accurately  using 
the  compound  Poisson  noise  of  polyenergetic  x-ray  beam 
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(c)  (d) 


(e)  (f) 

Fig.  7.  CBCT  of  the  anthropomoiphic  head  phantom:  (a)  analytical  FDK 
reconstructed  image  from  projection  data  acquired  using  a  low-dose  proto¬ 
col  (10  mA/ 10  ms)  and  (b)  a  high-dose  protocol  (80  mA/ 12  ms);  (c) 
PWLS  iterative  image  reconstruction  with  an  isotropic  quadratic  penalty 
from  projection  data  acquired  using  a  low-dose  protocol  and  (d)  with  a 
proposed  anisotropic  penalty,  (e)  PWLS  iterative  image  reconstruction  with 
the  Huber  penalty;  (f)  analytical  FDK  reconstructed  image  after  low-dose 
projection  data  processed  by  the  PWLS  criterion  (Ref.  34)  with  smoothing 
parameter  0.09.  Display  window;  [0,0.02]  mm“'. 


plus  Gaussian  electronic  noise. The  performance  of  itera¬ 
tive  reconstruction  algorithms  for  x-ray  CT  may  be  further 
improved  by  more  accurate  noise  modeling. 

The  penalty  reflects  the  prior  information  of  the  CT  im¬ 
ages.  In  this  work,  the  anisotropic  penalty  of  the  quadratic 
form  was  proposed  to  encourage  smoothness  among  neigh¬ 
boring  pixels  of  similar  intensities  but  discourage  the 
smoothness  if  a  large  difference  exists  between  neighboring 
pixels.  Thus,  edges  are  better  preserved  in  reconstructed  im¬ 
ages.  In  radiotherapy,  CT  images  of  the  same  patient  are 
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usually  available  before  treatment.  The  high-quality  planning 
CT  images  provide  strong  a  priori  information  of  the  patient 
and  it  may  be  used  to  improve  the  performance  of  iterative 
image  reconstruction  algorithms.  However,  interfractional 
variation  in  treatment  position  and  deformation  of  organs 
may  make  such  application  challenging.  Dedicated  registra- 
tion  algorithms  are  necessary  to  extract  information  from 
the  planning  CT  as  a  prior  constraint  for  iterative  reconstruc¬ 
tion. 

Based  on  the  same  noise  properties  of  projection  data,  the 
PWLS  objective  function  can  also  be  constructed  in  the  pro¬ 
jection  or  sinogram  domain  where  the  penalty  is  applied  be¬ 
tween  neighboring  projection  pixels. CT  images  can  then 
be  reconstructed  by  analytical  algorithms  such  as  FDK. 
Compared  with  fully  iterative  image  reconstruction  methods, 
the  strategy  of  projection  smoothing  followed  by  analytical 
image  reconstruction  is  advantageous  in  computational  effi¬ 
ciency  because  projection  and  backprojection  cycles  in  the 
iterative  image  reconstruction  algorithm  are  avoided.  Re- 
cently.  La  Riviere  and  Vargas  have  shown  potential  equiva¬ 
lence  between  the  image-domain  based  iterative  reconstruc¬ 
tion  methods  and  the  strategy  of  sinogram  restoration  plus 
analytical  filtered  backprojection  reconstruction.  The  studies 
performed  in  their  work  are  based  on  a  simple  isotropic 
quadratic  penalty.  It  will  be  interesting  to  perform  a  similar 
study  based  on  the  proposed  anisotropic  quadratic  penalty. 
The  edge-preserving  penalty  in  image  domain  may  have 
some  advantages  compared  with  the  same  penalty  used  in 
projection  domain  because  edges  are  better  defined  in  the 
image  domain  than  of  that  in  the  projection  domain.  In  this 
work,  we  also  included  the  results  obtained  using  the  strat¬ 
egy  presented  in  Ref.  34,  i.e.,  the  projection  image  is  pro¬ 
cessed  according  to  the  PWLS  criterion  before  the  analytical 
FDK  reconstruction.  Figures  3(f)  and  7(f)  show  the  results 
from  the  projection-domain  approach^"^  with  smoothing  pa¬ 
rameter  yS=0.09  for  the  CatPhan®  600  phantom  and  the  an¬ 
thropomorphic  head  phantom,  respectively.  It  can  be  ob¬ 
served  that  the  edges  in  the  image  reconstructed  by  FDK 
from  the  PWLS  processed  projection  image  are  blurred  com¬ 
pared  with  the  image  reconstructed  by  the  PWLS  using  the 
anisotropic  quadratic  penalty.  For  a  quantitative  comparison, 
we  calculated  the  MTF  and  noise  level  of  the  image  of 
CTP591  module.  The  MTF  curve  was  plotted  in  Fig.  5  and 
the  standard  deviation  around  the  uniform  region  was  0.74 
X  10“^.  The  MTF  curves  in  Fig.  5  show  that  PWLS  image 
reconstruction  using  the  anisotropic  quadratic  penalty  pro¬ 
duces  better  image  resolution  at  the  matched  noise  level. 
This  initial  comparison  study  indicates  that  the  edge¬ 
preserving  penalty  in  the  image  domain  produces  higher  im¬ 
age  resolution  than  the  same  penalty  applied  in  the  projection 
domain  because  better  edge  definition  is  in  the  image  domain 
than  the  projection  domain. 

In  this  work,  our  effort  was  focused  on  the  noise  suppres¬ 
sion  of  CBCT  using  the  iterative  reconstruction  algorithm. 
The  presented  iterative  reconstruction  algorithm  can  also  be 
used  to  improve  image  quality  of  the  4D-CBCT.^®’^’  In  4D- 
CBCT,  projection  views  for  a  specific  phase  are  usually  ir¬ 
regular  and  undersampled.  Direct  reconstruction  using  the 
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conventional  FDK  algorithm  from  phase-binned  projection 
data  may  lead  to  unacceptable  results  due  to  view  aliasing 
artifacts.  Several  strategies,  such  as  using  slow-rotating 
imager  and  interphase  motion  model,  ’  have  been  pro¬ 
posed  to  enhance  the  image  quality  of  4D-CBCT.  Iterative 
reconstruction  algorithms  incorporate  both  data  acquisition 
geometry  and  sampling  of  projection  views  into  the  projec¬ 
tion  matrix  automatically.  Consequently,  the  quality  of  the 
CBCT  image  so  obtained  is  generally  superior  over  that  re¬ 
constructed  using  an  analytical  method."^*' 

Although  iterative  reconstruction  algorithms  have  shown 
advantages  for  CT  imaging  in  terms  of  noise  suppression  and 
structure  preservation,  computational  time  could  be  a  chal¬ 
lenge  for  its  practical  use.  In  our  implementation,  we  com¬ 
puted  the  projection  matrix  A  before  iterative  reconstruction. 
The  projection  matrix  was  stored  as  a  hie  and  served  as  a 
lookup  table  during  iterations.  It  takes  about  15  min  to  hnish 
one  iteration  to  reconstruct  the  CBCT  images  of  a  size  350 
X  350  X  16  using  a  PC  with  3  GHz  CPU.  Nevertheless,  the 
reconstruction  can  be  sped  up  by  graphics  card 

acceleration*®’"^^  and  parallel  computation  using  PC  clusters"*^ 

18 

and  cell  broadband  engine. 


V.  CONCLUSION 

In  this  work,  we  presented  a  statistics-based  iterative  re¬ 
construction  algorithm  for  CBCT.  The  objective  function 
was  based  on  the  PWLS  criterion.  To  preserve  edges  in  the 
reconstructed  images,  an  anisotropic  quadratic  penalty  was 
proposed.  Noise  and  artifacts  in  low-dose  CBCT  are  greatly 
suppressed  using  the  presented  PWLS  reconstruction  algo¬ 
rithm.  Comparison  studies  with  reconstruction  based  on  the 
isotropic  penalty  have  clearly  shown  the  beneht  of  the  pro¬ 
posed  approach.  The  statistical  iterative  reconstruction  algo¬ 
rithm  signihcantly  improves  low-dose  CBCT  image  quality 
and  may  hnd  useful  clinical  applications  in  the  future. 
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APPENDIX:  IMAGE  RECONSTRUCTION 
ALGORITHM 

The  task  for  image  reconstruction  is  to  estimate  the  at¬ 
tenuation  coefficient  distribution  map  /x  from  the  projection 
data  p  by  minimizing  Eq.  (2).  In  this  study,  the  minimization 
was  performed  iteratively  using  the  Gauss-Seidel  update  al¬ 
gorithm,  similar  to  that  in  Ref.  25, 


Initialization: 


A=FDK{p} 

f=y-A(L 

X  =  diagjof  (/7,.)}= diag{exp(p,.)  /  A  J 


i-ajX-'fly,  Vj, 
exp  - 


niGNj 


For  each  iteration: 
begin 

For  each  pixel  j: 
begin 


*'  old . 


aj:=Sj+pi, 


m&N^jm 

;x-*f+.,yx°“*+;SX 


jlj  ■■=  max{0 ,  /x“"} 


Wjm-- 

end 


Wj„,  exp 


end 


(Al) 


The  iterations  can  be  stopped  by  setting  a  threshold  for 
the  change  of  objective  function  or  the  number  of  iterations. 
In  all  of  cases  presented  in  this  work,  we  stopped  the  com¬ 
putation  at  20  iterations  at  which  good  convergence  was  seen 
through  the  observation  of  the  reconstructed  image  at  each 
iteration. 
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Scatter  correction  is  crucial  to  the  quality  of  reconstructed  images  in  x-ray  cone-beam  computed 
tomography  (CBCT).  Most  of  existing  scatter  correction  methods  assume  smooth  scatter  distribu¬ 
tions.  The  high-frequency  scatter  noise  remains  in  the  projection  images  even  after  a  perfect  scatter 
correction.  In  this  paper,  using  a  clinical  CBCT  system  and  a  measurement-based  scatter  correction, 
the  authors  show  that  a  scatter  correction  alone  does  not  provide  satisfactory  image  quality  and  the 
loss  of  the  contrast-to-noise  ratio  (CNR)  of  the  scatter  corrected  image  may  overwrite  the  benefit  of 
scatter  removal.  To  circumvent  the  problem  and  truly  gain  from  scatter  correction,  an  effective 
scatter  noise  suppression  method  must  be  in  place.  They  analyze  the  noise  properties  in  the  pro¬ 
jections  after  scatter  correction  and  propose  to  use  a  penalized  weighted  least-squares  (PWLS) 
algorithm  to  reduce  the  noise  in  the  reconstructed  images.  Experimental  results  on  an  evaluation 
phantom  (Catphan©600)  show  that  the  proposed  algorithm  further  reduces  the  reconstruction  error 
in  a  scatter  corrected  image  from  10.6%  to  1.7%  and  increases  the  CNR  by  a  factor  of  3.6. 
Significant  image  quality  improvement  is  also  shown  in  the  results  on  an  anthropomorphic  phan¬ 
tom,  in  which  the  global  noise  level  is  reduced  and  the  local  streaking  artifacts  around  bones  are 
suppressed.  ©  2009  American  Association  of  Physicists  in  Medicine.  [DOI:  10.1118/1.3063001] 


I.  INTRODUCTION 

In  x-ray  computed  tomography  (CT),  scatter  causes  severe 

1  2 

distortions  and  contrast  loss  in  the  reconstructed  images.  ’ 
Scatter  magnitude  increases  as  the  x-ray  illuminated  volume 
size  increases.  In  an  x-ray  system  with  a  large  area  detector, 
such  as  a  cone-beam  CT  (CBCT)  system,  the  scatter-to- 
primary  ratio  (SPR)  can  be  as  high  as  8  in  certain  areas  of  the 
projection  images."^’^ 

CBCT  is  being  commonly  used  in  many  applications  for 
its  large  volume  coverage.  However,  the  high  SPR  severely 
deteriorates  the  quality  of  CBCT  image  and  hampers  its 
clinical  usage.  Many  scatter  correction  methods  have  been 
proposed  in  the  literature,  and  research  in  this  field  is  still 
very  active.^’®^^  There  are  two  major  types  of  scatter  removal 
techniques.  The  first  type  performs  scatter  suppression  dur¬ 
ing  the  acquisition  of  projection  data,  based  on  the  difference 
between  the  incident  angles  of  primary  photons  and  scatter 
photons.  Typical  examples  include  the  antiscatter  grid 
method  and  the  air  gap  method. Although  instant  scatter 
suppression  is  achieved  using  these  methods,  their  efficacy  is 
usually  limited.  Siewerdsen  et  al,  for  example,  showed  that 
an  antiscatter  grid  was  effective  only  in  improving  the 
contrast-to-noise  ratio  (CNR)  of  low  resolution  CT  images.^* 
Kyriakou  et  al.  also  reported  that  if  an  antiscatter  grid  is  used 
and  the  scatter  is  high,  the  imaging  dose  need  to  be  increased 
significantly  to  compensate  for  the  primary  loss  due  to  the 
insertion  of  the  grid.  An  improved  scatter  correction  capa¬ 
bility  has  been  demonstrated  using  the  scatter  correction 
methods  in  the  second  category,  where  the  scatter  is  cor¬ 
rected  using  postprocessing  techniques  on  the  scatter  con¬ 
taminated  projection  images. Due  to  the  randomness  of 
scattering  events,  however,  these  methods  implicitly  or  ex¬ 
plicitly  assume  smooth  scatter  distributions  and  tacitly  ignore 
the  existence  of  high-frequency  scatter  noise.  The  scatter 


noise  is  therefore  left  in  the  image  after  scatter  correction, 
resulting  a  degradation  of  CNRs  in  the  reconstructed 
images.^®  So  far,  the  noise  due  to  scatter  in  the  x-ray  projec¬ 
tion  image  is  generally  considered  as  a  low-dose  imaging 
problem,  and  little  attention  has  been  paid  in  the  literature  to 
reduce  the  noise  due  to  scatter  correction.  As  the  scatter  cor¬ 
rection  techniques  become  more  successful,  this  issue  is  be¬ 
coming  increasingly  important. 

The  purpose  of  this  paper  is  to  investigate  the  role  of 
high-frequency  noise  in  the  CB  projection  image  after  scatter 
correction  and  to  provide  a  practical  solution  to  noise  sup¬ 
pression  in  scatter  corrected  reconstructed  images.  We  first 
study  the  noise  property  of  the  projection  images  after  scatter 
correction  according  to  Poisson  statistics.  A  penalized 
weighted  least-squares  (PWLS)  algorithm  is  then  applied  to 
effectively  suppress  the  image  noise. The  algorithm  is 
evaluated  using  phantom  experiments  on  a  clinical  CBCT 
system. 

II.  METHOD 

II.A.  Noise  in  x-ray  projections  with  and  without 
scatter  correction 

There  are  two  major  types  of  noise  in  x-ray  projection 
images.  One  is  the  image  independent  noise  due  to  the  elec¬ 
trical  and  roundoff  error,  which  can  be  considered  as  Gauss¬ 
ian  noise;  the  other  is  the  image  dependent  noise  due  to  the 
statistical  fluctuation  of  the  x-ray  photons  that  exit  an  object, 
which  can  be  considered  as  Poisson  noise.  We  assume  the 
noise  of  the  first  type  is  small  and  consider  only  the  Poisson 
noise  here. 

Denote  variables  s  and  p  as  the  detected  scatter  and  pri¬ 
mary  signals  with  mean  values  S  and  P,  respectively,  s  and  p 
are  Poisson  distributed,  i.e.,  i~Poisson(5)  and  p 
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Fig.  1.  Reconstructed  images  of  the  Catphan©600  phantom  with  an  oval  body  annulus.  Display  window:  [-500  500]  HU.  The  mean  and  standard  deviation 
(std)  inside  the  white  squai'es  in  the  images  are  measured  as  mean+std  HU.  (a)  No  scatter  correction  and  no  noise  suppression;  CT  number  in  the  selected  ROI 
(white  squai'e):  -117  ±51  HU.  (b)  Scatter  coiTection  without  noise  suppression;  CT  number  in  the  ROI:  17  ±191  HU.  (c)  Scatter  correction  using  the 
proposed  noise  suppression  algorithm,  yS=0.0009;  CT  number  in  the  ROI:  18  ±  52  HU.  (d)  Scatter  correction  using  the  proposed  noise  suppression  algorithm, 
/5=0.0001;  CT  number  in  the  ROI:  16±  108  HU.  (e)  Scatter  correction  with  noise  suppression  using  the  standard  Hamming  filter;  CT  number  in  the  ROI: 
18  ±  152  HU.  (f)  Scatter  correction  using  the  proposed  noise  suppression  algorithm,  assuming  i'^=0;  CT  number  in  the  ROI:  18  ±  151  HU. 
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S  +  P 


(n,  +  n„). 
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var(«nJ 
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{S  +  Pf 
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var(«^  +  n„) 


{S  +  Pf 
1 


iS  +  P) 


S  +  P 


(5) 


When  an  effective  scatter  correction  algorithm  is  applied, 
the  scatter  mean  value  S  is  removed,  while  the  high- 
frequency  noise  is  left  in  the  corrected  image.  Similarly, 
the  scatter  corrected  line  integral  image  can  be  written  as 


Fig.  2.  Measured  MTFs  using  the  Catphan©600  phantom.  The  algorithm 
parameters  are  tuned  such  that  the  last  three  curves  as  shown  in  the  legend 
match. 


~  Poisson(P).  From  the  fact  that  the  point  spread  function  of 
the  scatter  is  broad  and  smooth,^®  we  assume  s  is  indepen¬ 
dent  of  p,  and  therefore,  (i-tp)  ~  Poisson(5'H-P).  Denote 
and  Hp  as  the  statistical  noises  of  the  scatter  and  primary, 
respectively,  and  one  obtains 


q,  =  Inl 


:ln 


s  +  p-  S 
Iq 

n,  +  np  +  P 


■  ln(u,  +  tip  +  P)  +  ln(/o) 

■  ln(P)  H-  ln(/o)  -^{ns  +  tip)  +  Oi{n,  +  Hp)^) .  (6) 


The  associated  noise  and  its  variance  is  approximated  as 


p  =  P  +  np. 

Due  to  Poisson  statistics. 


var{n^)  =  S, 


vaT{np)  =  P. 


(2) 


Denoting  Iq  as  the  incident  photon  intensity  and  q  as  the 
line  integral  image  without  scatter  correction,  we  have 


g'  =  In 

\s  +  p  / 

=  -ln(sH-p)H-ln(/o) 


1 


=  -  ln(5  +  P)  +  ln(/o)  -  - — -{n,  +  tip)  +  Oi{n,  +  Hp)^) . 

O  1 


(3) 


The  last  step  uses  Taylor’s  expansion  at  (S+P). 

Since  Iq  is  typically  very  large,  the  noise  associated  with 
log(/o)  can  be  ignored.  Assume  (n^+Hp)  is  small  and  ignore 
the  high  order  term,  then  the  noise  of  q,  and  its  variance 
are  approximated  as 


^  -  -itls  +  tlp) 


\ax{n^)  =  ^var(n^  H-  tip) 


JiiS  +  P) 


S  +  P 


S 

1  H-  - 
p 


(7) 


Therefore,  the  ratio  of  the  variances  of  and  is 


var(nc)  ^  L 
var(«„c)  \  pI 


(8) 


Equation  (8)  shows  that  the  noise  is  magnified  after  a 
scatter  correction  is  applied  on  the  projection  image.  Since 
the  noise-free  scatter  signal  S  is  spatially  much  smoother 
than  the  noise-free  primary  signal  P,  the  spatial  distribution 
of  SPR  (S/ P)  is  typically  nonuniform.  In  a  projection  image 
on  a  human  chest,  for  example,  the  SPR  value  is  very  low 
around  the  object  boundary  (less  than  0.2)  and  can  be  larger 
than  8  in  areas  of  small  primary  signals  (e.g.,  behind  bones). 
As  a  result,  the  noise  variance  magnification  across  the  pro¬ 
jection  image  ranges  from  below  1.5  to  more  than  80.  As  will 
be  shown  in  Sec.  Ill,  the  nonuniform  noise  magnification  not 
only  increases  the  noise  level  in  the  reconstructed  image  but 
also  results  in  local  streaking  artifacts. 
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(d) 

Fig.  3.  Difference  images  after  applying  noise  suppression  on  Fig.  1(b)  using  different  algorithms.  Display  window:  [-150  150]  HU.  (a)  Scatter  correction 
using  the  proposed  noise  suppression  algorithm,  /3=0.0009.  [Fig.  l(c)-Fig.  1(b)].  (b)  Scatter  coiTection  using  the  proposed  noise  suppression  algorithm,  /3 
=0.0001.  [Fig.  l(d)-Fig.  1(b)].  (c)  Scatter  correction  with  noise  suppression  using  the  standard  Hamming  filter  [Fig.  l(e)-Fig.  1(b)].  (d)  Scatter  correction 
using  the  proposed  noise  suppression  algorithm,  assuming  5^=0.  [Fig.  l(f)-Fig.  1(b)]. 


II.B.  The  noise  suppression  aigorithm 

Based  on  the  noise  property  of  the  scatter  corrected  image 
as  derived  in  Eq.  (7),  we  modify  and  implement  a  previously 
developed  penalized  weighted  least-squares  (PWLS) 
method  to  suppress  the  image  noise.  The  PWLS  method  is 
a  statistics-based  algorithm  that  aims  to  estimate  the  ideal 
linear  integrals  by  minimizing  the  PWLS  objective  function. 
The  PWLS  objective  function  models  the  first  and  second 
moments  of  the  projection  data.  Mathematically,  the  PWLS 
cost  function  is  written  as 

<Piqc)  =  iqc  -  -  Re)  +  I^Riqc)-  (9) 

The  first  term  in  Eq.  (9)  is  a  weighted  least-squares  criterion. 
The  variable  is  the  vector  of  the  scatter  corrected  line 
integral  data  as  shown  in  Eq.  (6);  the  variable  is  the  vector 
of  noise  suppressed  line  integral  data  to  be  estimated.  The 
symbol  T  denotes  the  transpose  operator.  The  matrix  S  is 
diagonal  matrix  and  its  /th  element  is  the  variance  of  q^  at 


detector  pixel  i.  Erom  Eq.  (7),  we  approximate  the  variance 
of  q^  based  on  the  following  equations: 

p,„=^S  +  P 


variq^)  =  var{n,) 


1  /  5 

« -  l-H- 

S  +  P\  P 
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Fig.  4.  Reconstruction  of  the  Catphan©600  phantom  using  the  proposed 
algorithm  with  an  accurate  scatter  estimate  in  the  scatter  con'ection  but  an 
inaccurate  scatter  estimate  in  the  noise  suppression,  (a)  Noise  suppression 
using  the  proposed  algorithm,  assuming  scatter  to  be  a  constant  fraction  of 
the  smallest  intensity  in  each  projection;  CT  number  in  the  ROI:  15  ±  183 
HU.  Display  window:  [-500  500]  HU.  (b)  Difference  image  [Fig.  4(a)-Fig. 
1(b)].  Display  window:  [-150  150]  HU. 


where  is  the  measured  projection  data  before  scatter  cor¬ 
rection  and  logarithm  operation;  is  the  scatter  estimate 
obtained  from  scatter  correction  algorithms.  With  the  non- 
uniform  weights  included  in  the  least-squares  criterion,  the 
data  with  lower  signal-to-noise  ratios  (SNR)  contribute  less 
for  estimation  of  the  noise  suppressed  line  integrals. 

The  second  term  in  Eq.  (9)  is  a  smoothness  penalty  or  a 
prior  constraint,  where  ^  is  the  smoothing  parameter  which 
controls  the  degree  of  agreement  between  the  estimated  and 
the  measured  data.  A  quadratic  function  as  used  in  the  pre- 
vious  publications  is  adopted  here, 

R{q)  =  '2,'2,Wi,iqi-q„f,  (12) 

i  n 

where  n  represents  four  nearest  neighbors  around  pixel  i  and 
Wi„  is  the  weight  for  neighbor  n.  To  preserve  the  edge  infor¬ 
mation,  anisotropic  weights  are  used  as  w„,.  The  weight  of 


the  neighbor  is  determined  by  the  magnitude  of  difference 
between  neighbors  and  the  pixel  of  concern. 


:  exp 


qi  -  qn 


(13) 


where  S  is  another  user-defined  parameter  which  controls  the 
strength  of  edge  preservation. 

In  our  implementations,  the  parameter  S  was  chosen  as 
such  a  value  that  90%  of  the  sinogram  pixels  in  the  projec¬ 
tion  images  to  be  processed  has  a  gradient  magnitude  smaller 
than  SP'  The  cost  function  (9)  is  minimized  efficiently  using 
an  iterative  Gauss-Seidel  updating  strategy.  If  standard 
CBCT  imaging  parameters  are  used,  the  algorithm  converges 
to  an  optimal  solution  typically  after  about  20  iterations.  On 
a  3.0  GHz  PC,  the  process  takes  about  3  s  on  each  projection 
image. 


II.C.  Scatter  estimation 

In  this  work,  we  use  a  measurement-based  method  for 
scatter  estimation.  In  each  experiment,  a  lead  strip  was  in¬ 
serted  between  the  x-ray  source  and  the  object,  resulting  in  a 
horizontal  strip  shadow  on  the  detector  with  a  width  of  ap¬ 
proximately  2  cm.  Since  the  lead  strip  attenuates  almost  all 
the  incident  photons,  the  signal  detected  inside  the  strip 
shadow  only  contains  scatter.  Assuming  that  the  scatter  dis¬ 
tributions  with  and  without  the  insertion  of  the  lead  strip  are 
approximately  the  same,  we  consider  the  measured  scatter 
signals  using  a  lead  strip  as  the  scatter  signals  in  a  conven¬ 
tional  scan  at  the  same  location.  To  get  rid  of  the  noise  in  the 
measurement,  the  detected  signals  inside  the  2  cm  strip 
shadow  were  first  averaged  in  the  longitudinal  direction  and 
then  smoothed  laterally  using  a  moving  average  filter.  The 
resultant  signal  was  used  as  the  scatter  estimate  in  a  follow¬ 
ing  conventional  scan. 


II. D.  Evaluation 

The  proposed  method  is  evaluated  using  experiments  on  a 
Varian  Acuity  CT  simulator  (Varian  Medical  Systems,  Palo 
Alto,  CA).  This  CBCT  imaging  system  is  commonly  used  in 
radiation  therapy  for  verifying  the  patient  geometry.  A  stan¬ 
dard  protocol  as  used  in  clinic  was  used  in  the  experiments. 
The  x-ray  tube  operated  at  125  kVp  voltage  and  80  mA  with 
the  pulse  width  at  each  projection  angle  of  25  ms.  Data  of  a 
360°  scan  consist  of  about  680  projections  with  an  angle 
interval  of  about  0.5°.  The  dimension  of  each  acquired  pro¬ 
jection  image  was  397.3  X  298.0  mm^,  containing  1024 
X  768  pixels.  The  source-to-axis  distance  (SAD)  was 
1000  mm  and  the  source-to-imager  distance  (SID)  was 
1500  mm.  To  increase  the  field  of  view  (FOV),  a  half-fan 
mode  was  used,  with  the  flat-panel  detector  shifted  by  ap¬ 
proximately  160  mm. 

Two  phantoms  were  used  in  the  experiments.  The  first 
was  an  evaluation  phantom  which  consisted  of  a  Cat- 
phan©600  phantom  in  the  middle  with  a  radius  of  10  cm  and 
an  oval  body  annulus  in  the  periphery  to  expand  the  phantom 
to  an  elliptical  cylinder  with  a  major  axis  of  38  cm  and  a 
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Fig.  5.  Reconstructed  images  of  the  Catphan©600  phantom  with  an  oval  body  annulus.  Display  window:  [-500  1700]  HU.  A  different  slice  from  Fig.  1  is 
shown  to  investigate  the  resolution  performance.  A  zoom-in  image  of  the  line  pairs  inside  the  dashed  rectangle  is  shown  at  the  lower  left  corner  in  each  image. 
The  dashed  arcs  indicate  the  sets  of  line  pairs,  of  which  ID  profiles  passing  through  the  center  are  compared  in  Fig.  6.  (a)  No  scatter  comection  and  no  noise 
suppression,  (b)  Scatter  correction  without  noise  suppression,  (c)  Scatter  correction  using  the  proposed  noise  suppression  algorithm,  yS=0.0009.  (d)  Scatter 
correction  using  the  proposed  noise  suppression  algorithm,  yS=0.0001. 


minor  axis  of  30  cm.  The  uniform  oval  annulus  is  made  of 
the  same  material  as  that  of  the  CTP  486  uniformity  test 
module  inside  the  standard  Catphan©600  phantom,  and  it 
has  an  estimated  CT  number  of  —15  Hounsfield  Units  (HU). 
Due  to  the  increased  volume  size,  scatter  was  high  in  the 
scan  and  the  image  quality  was  much  worse  than  that  when 
only  the  Catphan©600  phantom  was  imaged.  No  bow-tie 
filter  was  used  in  the  experiments.  The  second  phantom  was 
an  anthropomorphic  chest  phantom.  A  bow-tie  filter  as  used 
in  clinic  was  placed  on  in  order  to  maintain  a  more  uniform 
photon  statistics  across  the  FOV.  The  reconstructed  images 
are  presented  in  HU,  i.e.,  with  a  CT  number  of  -1000  HU 
for  air  and  a  CT  number  of  0  HU  for  water-equivalent  ma¬ 
terials. 

Images  are  reconstructed  using  the  standard  FDK 
algorithm,  without  scatter  correction  and  noise  suppression, 
with  scatter  correction  but  without  noise  suppression,  and 
with  scatter  correction  and  the  proposed  noise  suppression 
using  different  p  values.  The  same  filters  are  used  in  all  the 


reconstructions.  To  demonstrate  the  advantage  of  the  pro¬ 
posed  algorithm  on  noise  suppression,  two  more  reconstruc¬ 
tions  of  the  Catphan©600  phantom  are  carried  out  with  the 
same  scatter  correction.  For  a  fair  comparison,  the  noise  sup¬ 
pression  parameters  in  these  reconstructions  are  adjusted 
such  that  the  corresponding  modulation  transfer  functions 
(MTF)  approximately  match,  indicating  a  similar  spatial  res¬ 
olution  performance.  The  first  reconstruction  uses  a  standard 
low-pass  Hamming  filter  in  the  filtering  step  in  the  FDK 
algorithm  to  suppress  the  projection  noise.  The  second  re¬ 
construction  assumes  zero  scatter  in  the  noise  suppression 
step  after  the  scatter  correction  and  directly  uses  the  previ¬ 
ously  developed  PWLS  algorithm  for  the  noise  suppression, 
i.e.,  assuming  s^=Q  in  Eq.  (11).  The  tungsten  carbide  inside 
the  Catphan©600  phantom  with  a  diameter  of  0.28  mm  is 
used  for  the  MTF  measurements. 

Side-by-side  comparisons  of  reconstructed  images  are 
provided  to  illustrate  the  performance  on  noise  reduction  in 
the  proposed  algorithm  as  well  as  its  effect  on  the  image 
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Fig.  6.  Comparison  of  ID  profiles  passing  through  the  center  of  the  line  7.  Catphan  image  without  the  oval  body  annulus  used  as  the  bench- 

pairs  indicated  by  the  dashed  arcs  shown  in  Fig.  5.  mark.  A  nari'ow  collimator  was  also  used  to  further  reduce  the  scatter  arti¬ 

facts.  Seven  ROIs  are  used  in  the  quantitative  analysis. 


resolution.  Difference  images  before  and  after  noise  suppres¬ 
sion  are  also  presented.  In  the  quantitative  analysis,  besides 
image  noise  levels,  image  contrasts  and  CNRs  in  selected 
regions  of  interest  (ROI)  with  relatively  uniform  recon¬ 
structed  values  are  used.  The  image  contrast  is  defined  as 

contrast  =  /r,.  -  ,  (14) 


where  /u-,.  is  the  mean  reconstructed  value  inside  the  ROI  and 
is  the  mean  reconstructed  value  in  the  surrounding  area. 
The  CNR  is  defined  as 


CNR  = 


(15) 


where  is  the  variance  inside  the  ROI  and  a\  is  the  vari¬ 
ance  in  the  surrounding  area. 


III.  RESULTS 

Figure  I  shows  the  reconstructed  images  of  the  Cat- 
phan©600  phantom  with  an  oval  annulus.  The  image  distor¬ 
tion  is  obvious  in  Fig.  1(a),  where  no  scatter  correction  is 
applied.  Most  of  the  shading  artifacts  due  to  scatter  are  elimi¬ 
nated  using  the  measurement-based  scatter  correction,  as 
shown  in  Fig.  1(b).  However,  the  noise  is  magnified  in  the 
scatter  corrected  image  and  the  image  quality  is  degraded. 
Figures  1(c)  and  1(d)  are  the  images  with  noise  suppressed 
using  the  proposed  method.  To  investigate  the  performance 
of  the  algorithm  on  noise  reduction  and  image  resolution, 
two  different  values  are  used.  The  first  value  is  relatively 
large,  and  generates  an  image  [Fig.  1  (c)]  with  approximately 
the  same  noise  level  in  the  selected  ROI  as  that  in  the  image 
without  scatter  correction  and  noise  suppression.  The  im¬ 
provement  of  image  quality  is  significant.  A  small  (3  value  is 
used  in  the  second  implementation  of  the  noise  suppression 
algorithm.  This  value  achieves  a  less  smoothed  image  as 
seen  in  Fig.  1(d).  Figure  1(e)  is  the  result  using  a  standard 
Hamming  filter  on  the  projection  images  for  noise  suppres¬ 


sion.  Figure  1(f)  is  the  result  using  the  previously  developed 
PWLS  algorithm.  Equivalently,  we  assume  ^^=0  in  Eq.  (11) 
of  the  proposed  algorithm.  The  MTEs  of  the  six  reconstruc¬ 
tions  are  measured  and  shown  in  Eig.  2.  Note  that  we  adjust 
the  noise  suppression  parameters  such  that  Eigs.  l(d)-l(f) 
have  approximately  the  same  MTE  curves.  With  similar  per¬ 
formances  on  image  resolution,  Hamming  filtering  and  the 
previously  developed  PWLS  method  reduce  the  noise  stan¬ 
dard  deviation  in  the  selected  ROI  from  191  HU  to  approxi¬ 
mately  150  HU,  and  the  proposed  algorithm  is  able  to  further 
reduce  the  noise  standard  deviation  down  to  108  HU.  An 
improved  image  quality  using  the  proposed  algorithm  is  also 
evident  in  the  comparison  shown  in  Pig.  1.  Pigure  3  shows 
the  difference  images  of  the  same  slice  as  in  Pig.  1  before 
and  after  noise  suppression  using  different  algorithms  as  de¬ 
scribed  above. 

Pigure  1(f)  indicates  that  it  is  important  to  include  the 
scatter  estimate  in  the  proposed  noise  suppression  algorithm. 
Although  scatter  estimation/correction  algorithms  are  now 
becoming  more  and  more  successful,  it  is  still  challenging  to 
achieve  an  accurate  scatter  estimate  in  many  applications. 
Another  reconstruction  is  carried  out  to  study  the  sensitivity 
of  the  proposed  algorithm  performance  with  respect  to  the 
accuracy  of  the  scatter  estimation.  When  scatter  correction  is 
inaccurate,  residual  scatter  artifacts  are  dominant  as  com¬ 
pared  to  the  image  noise.  Por  a  better  illustration  of  the  al¬ 
gorithm  performance  on  noise  suppression,  we  use  an  accu¬ 
rate  scatter  estimate  for  the  scatter  correction.  In  the  noise 
suppression  step  using  the  proposed  algorithm,  we  assume 
that  scatter  is  uniform  over  the  whole  field  and  it  equals  a 
constant  fraction  of  the  smallest  intensity  in  each  projection. 
The  algorithm  parameters  are  tuned  such  that  the  corre¬ 
sponding  MTP  curve  matches  those  of  Pigs.  l(d)-l(f).  The 
reconstructed  image  and  its  difference  image  as  compared  to 
the  result  without  noise  suppression  [Pig.  1  (b)]  are  shown  in 
Pig.  4.  The  inaccurate  scatter  estimation  not  only  signifi¬ 
cantly  compromises  the  noise  suppression  capability  of  the 
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Table  I.  Contrast  comparison  of  the  selected  ROIs  (HU).  The  values  in  parentheses  are  the  relative  errors  in 
percentage  (absolute  values)  with  respect  to  the  benchmark  image  shown  in  Fig.  7. 


ROI 

1 

2 

3 

4 

No  scatter  correction  [Fig.  1(a)] 

83  (66.1%) 

302  (65.0%) 

-421  (60.7%) 

-98  (64.2%) 

Scatter  correction  only  [Fig.  1(b)] 

283  (15.7%) 

998  (13.3%) 

-1131  (5.5%) 

-279  (2.2%) 

Proposed  algorthum,  yS=  0.0009  [Fig.  1(c)] 

261  (6.4%) 

873  (1.1%) 

-1069  (0.2%) 

-272  (0.5%) 

Proposed  algorithm,  yS=0.0001  [Fig.  1(d)] 

267  (9.3%) 

876  (1.5%) 

-1067  (0.5%) 

-270  (1.3%) 

Using  a  Hamming  filter  [Fig.  1(e)] 

284  (15.9%) 

958  (11.0% 

-1123  (4.8%) 

-282  (3.0%) 

Proposed  algorithm,  5^=0  [Fig.  1(f)] 

281  (14.9%) 

956  (10.8%) 

-1123  (4.8%) 

-279  (1.9%) 

Benchmark  (Fig.  7) 

245 

863 

-1072 

-273 

ROI 

5 

6 

7 

mean  eiTor 

No  scatter  correction  [Fig.  1(a)] 

-61  (67.6%) 

-37  (70.0%) 

-334  (67.9%) 

65.9% 

Scatter  correction  only  [Fig.  1(a)] 

-213  (12.6%) 

-96  (22.8%) 

-1058  (1.9%) 

10.6% 

Proposed  algorithm,  yS=0.0009  [Fig.  1(c)] 

-190  (0.7%) 

-121  (2.8%) 

-1037  (0.1%) 

1.7% 

Proposed  algorithm,  yS=0.0001  [Fig.  1(d)] 

-202  (6.7%) 

-102  (18.3%) 

-1010  (2.7%) 

5.8% 

Using  a  Hmming  filter  [Fig.  1(e)] 

-209  (10.4%) 

-116  (6.9%) 

-1075  (3.5%) 

7.9% 

Proposed  algorithm,  5^=0  [Fig.  1(f)] 

-209  (10.8%) 

-104  (16.9%) 

-1065  (2.5%) 

8.9% 

Benchmark  (Fig.  7) 

-189 

-125 

-1038 

proposed  algorithm  [the  noise  standard  deviation  in  the  se¬ 
lected  ROI  of  Fig.  4(a)  is  183  HU]  but  also  results  in  addi¬ 
tional  streaking  artifacts  which  are  evident  in  the  difference 
image  [Fig.  4(b)].  It  is  interesting  to  note  that  the  proposed 
algorithm  with  a  zero  scatter  estimate  [Fig.  1(f)]  performs 
better  than  that  with  a  uniform  scatter  estimate  [Fig.  4(a)]. 
The  main  reason  is  that  the  performance  of  the  proposed 
algorithm  is  determined  by  the  estimation  accuracy  of  the 
spatial  distribution  of  the  image  noise,  instead  of  the  total 
magnitude.  Using  a  constant  fraction  of  the  smallest  intensity 
in  each  projection  as  the  scatter  estimate  gives  a  better  esti¬ 
mate  on  the  total  image  noise  magnitude  than  using  zero 
scatter  estimates.  However,  it  results  in  a  spatial  distribution 
of  noise  variance  more  different  from  the  truth.  More  discus¬ 
sion  on  this  issue  will  follow  in  a  later  section. 

To  further  investigate  the  image  resolution  with  respect  to 
the  (5  values.  Fig.  5  shows  the  reconstructed  slices  with  the 
resolution  gauges.  Note  that  a  display  window  different  from 
that  in  Fig.  1  is  used  for  a  better  display  of  the  line  pairs. 
Zoom-in  images  of  the  line  pairs  from  resolvable  to  irresolv¬ 
able  are  provided  at  the  lower-left  comer  of  each  image. 
Figure  6  shows  ID  profiles  that  pass  through  the  centers  of 
the  line  pairs  indicated  by  dashed  white  arcs  in  Fig.  5.  With 


(3  values  properly  chosen,  the  proposed  method  is  able  to 
significantly  reduce  the  noise  of  the  reconstructed  image 
with  a  negligible  resolution  loss. 

For  a  quantitative  analysis  of  the  image  quality,  another 
scan  was  carried  out  on  the  Catphan©600  phantom  without 
the  oval  body  annulus.  A  narrow  collimator  was  used  to  fur¬ 
ther  reduce  the  scatter.  The  reconstructed  image,  as  shown  in 
Fig.  7,  is  considered  as  a  benchmark  image.  Seven  ROIs  as 
indicated  in  Fig.  7  are  used  to  calculate  the  contrasts,  and  the 
results  are  summarized  in  Table  I.  The  error  relative  to  the 
contrasts  of  the  benchmark  image  is  provided  in  parentheses. 
The  scatter  correction  greatly  reduces  the  reconstruction  er¬ 
rors.  The  average  reconstruction  error  is  reduced  from  65.9% 
to  10.6%  when  the  scatter  correction  is  applied.  Using  the 
proposed  noise  suppression  algorithm,  the  average  error  is 
further  reduced  to  5.8%  when  a  small  /3  is  used  and  1.7% 
when  a  large  is  used.  The  comparison  also  shows  that 
when  the  image  resolution  is  matched,  the  proposed  algo¬ 
rithm  is  superior  to  the  standard  method  using  a  Hamming 
hlter  and  the  previously  developed  PWLS  algorithm. 

Table  II  shows  the  CNRs  for  different  ROIs.  After  scatter 
correction,  although  the  reconstruction  error  is  greatly  re¬ 
duced  and  image  contrast  is  increased,  the  average  CNR  de- 


Table  II.  CNR  comparison  of  the  selected  ROIs. 


ROI 

1 

2 

3 

4 

5 

6 

7 

Average 

No  scatter  coiTection  [Fig.  1(a)] 

1.37 

4.57 

6.57 

1.68 

1.28 

0.74 

5.78 

3.14 

Scatter  correction  only  [Fig.  1(b)] 

1.29 

3.94 

5.58 

1.15 

0.99 

0.38 

4.54 

2.55 

Proposed  algorithm,  yS=0.0009  [Fig.  1(c)] 

4.12 

14.16 

17.14 

4.87 

3.28 

2.09 

18.47 

9.16 

Proposed  algorithm,  yS=0.0001  [Fig.  1(d)] 

1.86 

6.07 

8.37 

2.08 

1.57 

0.74 

8.02 

4.10 

Using  a  Hamming  filter  [Fig.  1(e)] 

1.60 

4.90 

6.86 

1.42 

1.21 

0.57 

5.81 

3.19 

Proposed  algorithm,  5^=0  [Fig.  1(f)] 

1.58 

4.97 

6.98 

1.43 

1.14 

0.48 

5.59 

3.17 
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(b)  (d) 


Fig.  8.  Reconstructed  images  of  the  anthropomorphic  phantom.  Display  window:  [-500  500]  HU.  The  mean  and  standard  deviation  (std)  inside  the  white 
squares  in  the  images  are  measured  as  mean ±  std  HU.  (a)  No  scatter  correction  and  no  noise  suppression;  CT  number  in  the  selected  ROl  (white  square): 
-291  ±43  HU.  (b)  Scatter  correction  without  noise  suppression;  CT  number  in  the  selected  ROl:  1  ±  146  HU.  (c)  Scatter  conection  using  the  proposed  noise 
suppression  algorithm,  yS=0.001;  CT  number  in  the  selected  ROl:  17 ±46  HU.  (d)  Scatter  conection  using  the  proposed  noise  suppression  algorithm,  /3 
=0.0001;  CT  number  in  the  selected  ROl:  18  ±89  HU. 


creases  from  3.14  to  2.55  due  to  the  image  noise  magnifica¬ 
tion.  Using  the  proposed  noise  suppression  algorithm,  the 
average  CNR  increases  to  4.10  for  a  small  fS  and  9.16  for  a 
large  fS. 

Figure  8  shows  the  reconstruction  results  of  the  anthropo¬ 
morphic  phantom.  Note  that  a  bowtie  filter  was  used  in  the 
scan,  which  made  the  SPR  more  uniform.  However,  along 
the  projection  lines  which  pass  through  highly  attenuating 
objects,  such  as  bones,  the  SPR  can  still  be  very  high  (around 
8)  due  to  the  extremely  low  primary  signals.  After  the  loga¬ 
rithm  operation,  the  noise  in  those  projection  data  boosts  up 
after  an  effective  scatter  correction  and  results  in  streaking 
artifacts  in  the  reconstruction,  which  is  obvious  in  Fig.  8(b). 
As  shown  in  Figs.  8(c)  and  8(d),  these  artifacts  are  effec¬ 
tively  suppressed  using  the  proposed  algorithm,  and  the  im¬ 
age  noise  is  also  reduced.  Figure  9  shows  the  difference  im¬ 


ages  before  and  after  the  proposed  noise  suppression.  Figure 
10  compares  ID  central  vertical  profiles  of  the  images  shown 
in  Fig.  8. 

IV.  DISCUSSION 

Majority  of  the  existing  scatter  correction  methods  use 
postprocessing  techniques  on  the  scatter  contaminated  pro¬ 
jection  images  to  improve  the  reconstruction  accuracy.  These 
methods  only  partially  solve  the  problem  caused  by  scatter 
and  leave  the  high  frequency  components  of  the  scatter  sig¬ 
nal  intact,  which  often  leads  to  degraded  CNRs  and  image 
quality.  In  this  work,  both  theoretical  analysis  and  physical 
experiments  show  the  effect  of  noise  magnification  in  the 
reconstructed  image  due  to  scatter  correction.  Experiments 
were  carried  out  on  a  clinical  CBCT  system  with  a  com¬ 
monly  used  protocol  on  phantoms  with  a  human  size.  The 
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(b) 


Fig.  9.  Difference  images  after  applying  noise  suppression  on  Fig.  8(b) 
using  different  algorithm  parameters.  Display  window:  [-150  150]  FIU.  (a) 
Scatter  conection  using  the  proposed  noise  suppression  algorithm,  yS 
=0.001  [Fig.  8(c)-  Fig.  8(b)].  (b)  Scatter  coiTection  using  the  proposed  noise 
suppression  algorithm,  yS=0.0001.  [Fig.  8(d)-Fig.  8(b)]. 


results  indicate  that  an  effective  scatter  correction  alone  does 
not  provide  satisfactory  images  in  CBCT  because  the  gain 
from  the  scatter  reduction  is  inevitably  accompanied  with 
overwhelming  noise-related  artifacts.  One  traditional  method 
to  deal  with  the  noise  problem  is  to  increase  the  dose.  How¬ 
ever,  the  results  on  both  the  evaluation  phantom  and  the  an¬ 
thropomorphic  phantom  show  that  the  noise  variance  in¬ 
creases  by  a  factor  of  more  than  10  when  scatter  correction  is 
applied.  Since  noise  variance  of  a  CT  image  is  roughly  in¬ 
versely  proportional  to  the  total  number  of  incident 
photons,^"^  the  dose  needs  to  be  increased  by  a  factor  of  more 
than  10  in  order  to  make  the  noise  level  of  the  scatter  cor¬ 
rected  image  comparable  to  that  without  scatter  correction. 
This  excessive  dose  increase  is  hardly  acceptable  in  clinical 
practice. 

An  improved  PWLS  algorithm  is  implemented  in  this 
work  to  suppress  the  noise  in  reconstructed  images.  The 
PWLS  objective  function  models  the  first-order  and  second- 


Fig.  10.  Comparison  of  ID  central  vertical  profiles  of  the  images  shown  in 
Fig.  8.  The  results  using  the  proposed  noise  suppression  algorithm  with 
different  /5  values  have  a  relatively  large  difference  only  in  the  noisy  areas. 


order  statistics  of  the  measurements  and  it  is  equivalent  to 
the  penalized  maximum  likelihood  (pML)  or  maximum  a 
posteriori  (MAP)  criterion  for  the  independent  Gaussian  dis¬ 
tributed  noise.  Indeed,  the  noise  of  CT  line  integrals  can  be 
well  approximated  by  Gaussian  noise  based  on  an  experi¬ 
mental  study.^®  Therefore,  the  minimization  of  the  PWLS 
objective  function  gives  an  optimal  solution  in  a  statistical 
sense.  In  the  previous  work,  we  have  shown  that  the  perfor¬ 
mance  of  the  PWLS  algorithm  is  better  than  those  low-pass 
filters  during  image  reconstruction  or  noise  reduction  filters 
based  on  local  statistics  of  measurements.  In  this  paper, 
both  results  on  the  evaluation  phantom  and  the  anthropomor¬ 
phic  phantom  show  effective  reduction  in  global  image  noise 
as  well  as  local  streaking  artifacts  around  high  attenuating 
objects,  such  as  bones.  While  a  measurement-based  scatter 
correction  method  is  used  to  demonstrate  the  nature  of  the 
noise  magnification  problem  incurring  in  the  process  of  scat¬ 
ter  removal,  the  developed  noise  suppression  algorithm  is 
expected  to  work  effectively  with  other  scatter  correction 
algorithms  as  long  as  they  can  provide  accurate  scatter  esti¬ 
mates. 

One  concern  about  a  noise  suppression  algorithm  is  the 
possible  resolution  loss  due  to  smoothing.  We  compare  the 
proposed  algorithm  with  the  standard  noise  suppression 
method  using  a  Hamming  filter  on  the  projection  images. 
The  comparison  shows  that  when  a  similar  image  resolution 
is  achieved,  the  proposed  algorithm  is  superior  on  noise  sup¬ 
pression.  Our  results  also  show  that  using  a  conservative 
smoothing  strategy,  significant  noise  reduction  is  still  achiev¬ 
able  with  negligible  resolution  loss.  The  choice  of  parameter 
yS  in  the  PWLS  algorithm  is  a  tradeoff  between  reconstruc¬ 
tion  accuracy  and  resolution  loss.  In  adaptive  radiation 
therapy  using  CBCT,  for  example,  the  CBCT  images  are 
used  for  dose  calculation  and  the  CT  number  accuracy  is 
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more  important  than  the  image  resolution.  In  such  an  appli¬ 
cation,  an  aggressive  noise  suppression  strategy  (a  large  /3) 
should  be  used. 

The  proposed  algorithm  is  based  on  the  previously  devel¬ 
oped  PWLS  algorithm.  To  show  the  signihcance  of  including 
the  scatter  estimate  Sg  in  the  the  noise  estimation  formula 
[Eq.  (11)],  we  compare  the  results  using  the  proposed  new 
algorithm  and  the  previously  developed  PWLS  algorithm. 
The  latter  is  equivalent  to  the  proposed  algorithm  with  a  zero 
scatter  estimate.  Similar  to  the  comparison  with  the  result 
using  Hamming  filtering,  the  proposed  algorithm  achieves 
lower  noise  when  the  same  image  resolution  is  maintained. 
Note  that  this  degraded  ability  of  noise  suppression  using  the 
previously  developed  PWLS  algorithm  is  not  due  to  the  un¬ 
derestimation  of  total  noise  magnitude  by  assuming  zero 
scatter  signals,  since  such  an  underestimation  can  always  be 
compensated  for  using  a  larger  fS  in  the  algorithm.  The  es¬ 
sence  of  the  PWLS  algorithms  is  to  equalize  the  noise  vari¬ 
ance  of  different  pixels  by  assigning  different  weights.  The 
presence  of  scatter  greatly  changes  the  spatial  distribution  of 
the  noise  variance.  Since  the  noise  variance  is  proportional  to 
(l+SPR)^  according  to  Eq.  (10),  and  SPR  spatially  varies 
from  values  close  to  zero  to  those  larger  than  8,  the  true 
spatial  distribution  of  the  noise  variance  is  quite  different 
from  that  if  the  scatter  is  assumed  to  be  zero.  The  misesti- 
mation  of  scatter  results  in  incorrect  estimation  of  noise  spa¬ 
tial  distribution,  and  the  corresponding  weighting  in  the 
PWLS  algorithm  is  not  able  to  equalize  the  noise  variance  of 
each  projection  pixel.  The  noise  suppression  ability  is  there¬ 
fore  significantly  compromised.  To  illustrate  the  importance 
of  accurate  scatter  estimation  for  a  superior  noise  suppres¬ 
sion  performance,  we  also  compare  the  result  using  the  pro¬ 
posed  algorithm  with  a  uniform  scatter  estimate.  Since  uni¬ 
form  scatter  estimation  results  in  an  estimate  of  SPR 
distribution  much  less  accurate  than  that  using  a  zero  scatter 
estimate,  the  noise  suppression  performance  is  also  much 
worse.  Linally,  we  want  to  emphasize  that  the  inaccuracy  in 
scatter  estimation  compromises  the  performance  of  noise 
suppression  using  the  proposed  method,  however,  it  does  not 
make  reconstruction  less  accurate  than  that  without  noise 
suppression.  We  can  always  tune  the  yS  value  based  on  the 
tradeoff  between  the  noise  suppression  capability  and  the 
image  artifacts  caused  by  inaccurate  scatter  estimation. 

V.  CONCLUSION 

Scatter  correction  methods  based  on  postprocessing  fall 
short  in  eliminating  the  high-frequency  scatter  noise  and  are 
incapable  of  providing  satisfactory  CBCT  image  quality.  Us¬ 
ing  a  clinical  CBCT  system  with  conventional  imaging  set¬ 
tings,  we  have  shown  that  after  scatter  correction,  the  noise 
variance  in  selected  ROIs  of  the  reconstructed  image  can  be 
increased  by  a  factor  of  more  than  10.  We  argue  that  a  scatter 
correction  algorithm  should  be  used  together  with  a  noise 
suppression  algorithm  to  achieve  a  satisfactory  image.  A 
PWLS  algorithm  is  proposed  for  the  noise  suppression.  The 
phantom  experiments  indicate  that  the  algorithm  is  able  to 
further  reduce  the  reconstruction  error  in  a  scatter  corrected 


image  without  noise  suppression  from  10.6%  to  1.7%,  and 
the  average  CNR  in  selected  ROIs  is  also  increased  by  a 
factor  of  3.6.  The  improvement  in  image  quality  is  critical  in 
many  clinical  applications  of  CBCT,  such  as  accurate  dose 
calculation  and  tumor  target  delineation  in  radiation  therapy 
planning. 
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NOISE  REDUCTION  IN  LOW-DOSE  X-RAY  FLUOROSCOPY  FOR  IMAGE-GUIDED 
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Purpose:  To  improve  the  quality  of  low-dose  X-ray  fluoroscopic  images  using  statistics-based  restoration  algorithm 
so  that  the  patient  flnoroscopy  can  be  performed  with  reduced  radiation  dose. 

Method  and  Materials:  Noise  in  the  low-dose  fluoroscopy  was  suppressed  by  temporal  and  spatial  filtering.  The 
temporal  correlation  among  neighboring  frames  was  considered  by  the  Karhunen-Loeve  (KL)  transform  (i.e., 
principal  component  analysis).  After  the  KL  transform,  the  selected  neighboring  frames  of  fluoroscopy  were 
decomposed  to  uncorrelated  and  ordered  principal  components.  For  each  KL  component,  a  penalized  weighted 
least-squares  (PWLS)  objective  function  was  constructed  to  restore  the  ideal  image.  The  penalty  was  chosen  as 
anisotropic  qnadratic,  and  the  penalty  parameter  in  each  KL  component  was  Inversely  proportional  to  its  corre¬ 
sponding  eigenvalue.  Smaller  KL  eigenvalue  is  associated  with  the  KL  component  of  lower  signal-to-noise  ratio 
(SNR),  and  a  larger  penalty  parameter  should  be  used  for  such  KL  component.  The  low-dose  fluoroscopic  Images 
were  acquired  using  a  Varian  Acuity  simulator.  A  quality  assurance  phantom  and  an  anthropomorphic  chest  phan¬ 
tom  were  nsed  to  evaluate  the  presented  algorithm. 

Results:  In  the  images  restored  by  the  proposed  KL  domain  PWLS  algorithm,  noise  is  greatly  suppressed,  whereas 
fine  structures  are  well  preserved.  Average  improvement  rate  of  SNR  is  75%  among  selected  regions  of  interest. 
Comparison  studies  with  traditional  techniques,  such  as  the  mean  and  median  filters,  show  that  the  proposed 
algorithm  is  advantageous  in  terms  of  structure  preservation. 

Conclusions:  The  proposed  noise  reduction  algorithm  can  significantly  Improve  the  qnality  of  low-dose  X-ray  fluo¬ 
roscopic  image  and  allows  for  dose  reduction  in  X-ray  fluoroscopy.  ©  2009  Elsevier  Inc. 

Low-dose  fluoroscopy.  Penalized  weighted  least-sqnares,  Karhunen-Loeve  (KL)  transform.  Noise  reduction. 
Anisotropic  penalty. 


INTRODUCTION 

X-ray  fluoroscopic  imaging  plays  an  important  role  in  image- 
guided  radiation  therapy  (IGRT).  There  is  growing  interest  in 
using  X-ray  fluoroscopy  for  the  management  of  organ  motion 
and  gating  in  radiotherapy  (1^).  Both  room-mounted  or¬ 
thogonal  X-ray  fluoroscopic  imaging  system  (2)  and  gan¬ 
try-mounted  fluoroscopy  (3)  have  been  developed  for 
tumor  tracking  and  gating  for  radiotherapy.  Recently, 
a  method  for  real-time  tracking  of  implanted  hducial  marker 
was  developed  using  combination  of  kV  and  MV  imaging 
(4).  The  X-ray  fluoroscopic  imaging  provides  additional  in¬ 
formation  of  tumor  and  patient  structure;  however,  it  is  not 
risk-free.  The  extra  exposure  to  X-rays  during  fluoroscopy 
may  lead  to  adverse  health  effects  in  patients.  Therefore, 
low-dose  fluoroscopic  imaging  is  desirable  in  clinics. 

Low-dose  fluoroscopy  can  be  achieved  using  a  lower  mA 
and  pulse  length  during  acquisition  images  (5).  With  a  lower 


mAs  acquisition  protocol,  the  image  quality  will  be  degraded 
because  of  excessive  photon  quantum  noise.  Several  spatial 
and  temporal  noise  reduction  algorithms  (6-11)  have  been 
studied  to  remove  signal-dependent  noise  in  the  low-dose 
fluoroscopic  images.  Temporal  filtering  techniques  incorpo¬ 
rate  information  from  neighboring  frames  and  can  potentially 
improve  the  performance  of  two-dimensional  (2D)  spatial  fil¬ 
tering.  However,  temporal  filters  could  introduce  motion 
artifacts  in  the  filtered  fluoroscopic  images  if  motion  in  the 
fluoroscopy  is  not  handled  properly.  One  way  to  avoid  mo¬ 
tion  artifacts  resulting  from  temporal  filters  is  to  estimate 
the  motion  information  from  the  fluoroscopic  image 
sequences  and  compensate  for  it  during  filtering  process. 
Nevertheless,  estimation  of  motion  fields  from  2D  image 
sequences  is  challenging  because  of  its  ill-posed  nature  and 
excessive  noise  (12-14).  In  this  work,  we  aimed  to  reduce 
noise  at  low  doses  through  improved  temporal  filtering  and 
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statistics-based  image  restoration.  Instead  of  estimating  ex¬ 
plicit  motion  fields  from  2D  image  sequences,  we  considered 
the  correlation  between  neighboring  frames  of  fluoroscopy 
by  using  the  Karhunen-Loeve  (KL)  transform.  The  KL  trans¬ 
form  has  been  proved  to  be  useful  in  image  restoration  (15) 
and  tomographic  image  reconstruction  (16,  17).  The  KL 
transform  considers  similarities  as  well  as  differences 
between  neighboring  frames.  For  each  KL  component,  the 
noise  was  reduced  according  to  the  penalized  weighted 
least-squares  (PWLS)  criterion.  The  presented  method  was 
evaluated  using  two  phantom  experimental  studies. 

METHODS  AND  MATERIALS 

In  this  section,  we  first  introduce  the  PWLS  criterion  for  spatial 
noise  filtering.  We  then  describe  the  KL  transform  for  temporal  fil¬ 
tering,  followed  by  the  summary  of  the  implementation  of  the  KL 
domain  PWLS  algorithm.  We  later  describe  the  data  acquisition  of 
two  experimental  studies. 


Fig.  1.  Linear  relationship  between  mean  and  variance  of  X-ray 
fluoroscopic  image. 


Spatial  filtering  using  the  PWLS  criterion 

Mathematically,  the  PWLS  objective  function  can  be  described 
by  the  following  equation: 


-1 

(p{u}  =  [y  -  uf  {y  -  u)  +  ^R{u) .  (1) 

The  first  term  in  Eq.  1  is  a  weighted  least-squares  (WLS)  criterion, 
where  y  is  the  vector  of  the  measured  data  and  u  is  the  vector  of  ideal 
fluoroscopy  to  be  estimated.  T  denotes  the  transpose  operator.  The 
matrix  ^  is  the  covariance  matrix  of  y  and  its  ;th  diagonal  element 
a^i  is  the  variance  of  measured  data  at  detector  bin  i;  it  determines 
the  contribution  of  each  measurement  to  the  objective  function  be¬ 
cause  it  plays  the  weighting  role  in  the  WLS  criterion.  The  variance 
of  each  measurement  can  be  estimated  from  the  mean-variance  rela¬ 
tionship  of  Fig.  1.  To  obtain  the  mean-variance  relationship  of  low- 
dose  fluoroscopic  images,  600  repeated  fluoroscopic  images  of  the 
anthropomorphic  chest  phantom  were  acquired  when  the  platform 
was  static.  Mean  and  variance  of  each  pixel  were  then  calculated 
from  repeated  measurements  and  plotted  in  Fig.  1 .  It  can  be  observed 
that  mean  and  variance  of  measured  data  are  linear,  which  reflects 
the  Poisson  noise  nature  of  X-ray  photons.  A  non-zero  intersect 
can  also  be  observed  and  is  caused  by  background  electronic  noise. 

The  second  term  in  Eq.  1  is  a  smoothness  penalty  or  a  priori  con¬ 
straint,  where  jS  is  the  smoothing  parameter  that  controls  the  degree 
of  agreement  between  the  estimated  and  the  measured  data.  A  com¬ 
monly  used  penalty  is  a  quadratic  term  (17,  18): 

R{p)  =  plRp.  =  X] 

j  me  Nj 

where  index  j  runs  over  all  elements  in  the  fluoroscopic  image,  and 
Nj  represents  the  set  of  neighbors  of  the  /-th  pixel.  The  parameter  kj„, 
indicates  the  relative  contribution  of  different  neighbors  and  is  usu¬ 
ally  set  to  1  for  first-order  neighbors  and  1  / y/2  for  second-order 
neighbors.  One  drawback  of  this  choice  is  that  it  considers  only  dis¬ 
tance  information  of  neighbors  to  regularize  the  solution  and  may 
lead  to  oversmoothing  in  the  restored  image,  especially  around 
sharp  edges.  To  avoid  oversmoothing  around  edges,  we  chose  an 
edge-preserving  penalty  (19)  by  incorporating  the  gradient  informa¬ 
tion  into  the  parameter  kj„,.  The  weight  can  be  written  as: 


Km  =  Km  exp 


(3) 


This  new  parameter  has  a  form  similar  to  the  anisotropic  diffusion 
filter  (20),  in  which  the  gradient  and  the  parameter  5  determines  the 
strength  of  diffusion  during  each  iteration.  The  parameter  5  provides 
local  control  of  smoothness  and  can  be  chosen  as  a  value  such  that 
90%  of  the  pixels  in  the  image  to  be  processed  have  a  gradient  mag¬ 
nitude  smaller  than  5.  The  parameter  is  small  if  the  gradient  be¬ 
tween  the  neighbor  and  the  concerned  pixel  is  large.  The  neighbors 
with  large  gradient  usually  occur  at  the  edges,  and  equivalence  be¬ 
tween  such  neighbors  is  discouraged  by  introducing  Eq.  3.  There¬ 
fore,  the  edges  will  be  better  preserved  in  the  de-noised  image. 

The  minimization  of  the  PWLS  objective  function  (1)  can  be  per-  q2 
formed  efficiently  using  Gauss-Seidel  updating  strategy  (21): 

y>  +  M[  E  E 

(n+l)  _  _ \m^N; _ m^Nf 

1  +  /SVV;  E  Km 

meNi 

where  index  n  represents  iterative  number,  A,’  denotes  the  two  near¬ 
est  neighbors  of  i  in  which  the  index  is  smaller  than  ;,  denotes 
those  two  nearest  neighbors  of  i  in  which  the  index  is  larger  than 
i,  and  A,  denotes  these  four  nearest  neighbors  of  pixel  i  in  the  fluo¬ 
roscopic  image.  The  initial  of  u  is  given  by  the  measured  data  y. 


Temporal  filtering  using  the  KL  transform 

In  fluoroscopy,  different  image  frames  are  highly  correlated  for 
the  same  patient,  especially  among  neighboring  frames.  One  way 
to  use  the  correlation  among  neighboring  frames  is  the  KL  transform 
(KLT).  Following  KLT,  the  selected  neighboring  fluoroscopy 
frames  were  decomposed  to  uncorrelated  and  ordered  principal 
components.  Each  KL  component  was  associated  with  an  eigen¬ 
value,  which  can  be  effectively  used  during  spatial  filtering  of  fluo¬ 
roscopic  imaging.  In  this  work,  two  frames  of  the  image  preceding 
the  frame  under  consideration  were  chosen  to  perform  KLT.  This 
selection  of  neighboring  frames  is  based  on  idea  that  only  previously 
obtained  frames  are  available  to  process  the  chosen  frame  in  real 
time.  The  KL  transform  is  defined  as: 
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(a) 


(b) 


(c)  (d) 


Fig.  2.  Fluoroscopic  image  of  the  quality  assurance  phantom:  (a)  acquired  with  X-ray  tube  current  10  mA  and  pulse  length  2 
ms;  (b)  image  (a)  processed  by  the  Karhunen-Loeve  (KL)  domain  penalized  weighted  least-squares  (PWLS)  noise  reduction 
algorithm;  (c)  acquired  with  X-ray  tube  current  10  mA  and  pulse  length  10  ms.  (d)  Image  in  panel  a  processed  with  the  KL 
domain  PWLS  using  an  isotropic  quadratic  term  as  a  penalty.  In  panel  a,  b  1,  b2,  b3 ,  and  b4  indicate  the  background  region  for 
the  calculation  of  signal-to-noise  ratios.  Blurred  edges  can  be  observed  in  panel  d.  ROI  =  region  of  interest. 


L=Ay,  (5) 

where  y  is  vector  of  selected  frames  of  fluoroscopy  and  is  arranged  in 
a  such  way  that  each  row  comprised  one  frame  of  the  fluoroscopic 
image,  y  is  the  KL-transformed  fluoroscopic  image,  and  A  is  the 
KLT  matrix  that  can  be  calculated  from  y.  From  the  selected  frames 
of  fluoroscopy  y,  the  elements  of  the  covariance  matrix  K,  can  be  cal¬ 
culated  by: 

1  ® 

[k<]u=  -  yi)  ’ 

/=i 

where  is  the  datum  of  fluoroscopy  y  at  detector  bin  i  of  k-th  frame, 
y,-  /  is  the  datum  of  y  at  detector  bin  i  of  l-th  frame,  yk  yj  are  the 
mean  of  fluoroscopic  image  of  k-th  frame  and  /-th  frame  respectively. 


B  is  the  number  of  detector  bins  for  each  fluoroscopy  and  index  k  or  / 
runs  over  the  selected  nearby  frames.  From  the  covariance  matrix  K„ 
the  KLT  matrix  A  can  then  be  calculated  based  on 

K,A'  =  A'D.  (7) 

In  Eq.  1,  D  =  diag{di}]^^,  where  dj  is  the  /-th  eigenvalue  of  K,. 
After  the  KLT,  the  covariance  K,  of  y  will  be: 

k,  =  AK,A' =  diag{di}]^^.  (8) 

Equation  8  implies  that  the  covariance  matrix  of  the  KL-trans- 
formed  data  is  diagonal,  that  is,  the  covariance  of  the  signal  between 
different  frames  after  the  KLT  will  be  zero.  Therefore,  the  data  sig¬ 
nals  of  different  KL  components  are  no  longer  correlated  so  that 
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Fig.  3.  Profiles  across  the  bars  in  Fig.  2  (indicated  by  a  dashed  line 
in  Fig.  2a).  KL  =  Karhunen-Loeve;  PWLS  =  penalized  weighted 
least-squares. 


PWLS  objective  function  can  be  constructed  for  each  KL  principal 
component  separately.  For  each  principal  component  in  the  KL  do¬ 
main,  the  PWLS  cost  function  can  be  expressed  as  (16-18): 


=  {yi  -  '(y/-“/)  +  {P/di)R{ui), 


(9) 


where  yi  and  Ui  are  the  /-th  KL  components  of  y  and  u  respectively, 
and  is  the  diagonal  weighting  matrix  of  in  the  KL  domain  which 
can  be  estimated  by: 


(10) 


where  Qi  =  diag{aj^.  is  the  variance  matrix  of  the  projection  at  bin 

/,  is  the  variance  of  y,-  and  (p/  is  the  /-th  KL  basis  vector.  R{ui)  is  the 
edge-preserving  penalty  term  which  can  be  defined  as: 


R{u,)  =  ufid,  = 


^m,l ) 


(11) 


Equation  9  shows  that  the  smoothing  parameter  in  the  KL  domain 
shall  be  chosen  as  (^/4/),  and  the  PWLS  criterion  for  each  KL  com¬ 
ponent  becomes  adaptive  to  its  corresponding  eigenvalue.  This 
choice  is  favorable  because  the  smoothing  parameter  varies 
adaptively  according  to  the  signal-to-noise  ratio  (SNR)  of  that  com¬ 
ponent.  A  smaller  KL  eigenvalue  is  usually  associated  with  a  compo¬ 
nent  having  a  lower  SNR  (16)  and,  therefore,  a  larger  smoothing 
parameter  is  used  to  penalize  this  noisier  data  component. 

In  summary,  implementation  of  the  presented  KL  domain  PWLS 
noise  reduction  approach  for  the  low-dose  fluoroscopy  can  be 
described  as  follows: 

1.  For  a  /-th  frame  fluoroscopic  image,  the  (/-l)-th  and  the  (/'-2)-th 
frame  of  the  images  are  selected  for  the  KL  transform. 

2.  Compute  the  covariance  matrix  kf  from  the  selected  frames  of  the 
projections. 

3.  Calculate  the  KL  transform  matrix  A  according  to  Eq.  7. 

4.  Apply  the  KL  transform  on  the  selected  frames  of  fluoroscopic 
images. 

5.  Perform  PWLS  minimization  on  each  KL  component. 

6.  Apply  inverse  KL  transform  on  the  processed  KL  components 
for  the  estimate  of  the  ideal  fluoroscopy  at  the  chosen  j-th  frame. 


Table  1.  SNRs  of  different  ROIs  in  figure  2. 


ROIl 

ROE 

ROE 

ROM 

10  ms 

9.7 

10.5 

14.3 

9.9 

2  ms 

5.8 

3.8 

6.2 

5.3 

KL  domain  PWLS  2  ms 
,3=  200 

8.3 

7.9 

12.8 

7.5 

Abbreviations:  KL  =  Karhunen-Loeve;  PWLS  =  penalized 
weighted  least-squares;  ROI  =  region  of  interest. 


Data  acquisition 

The  X-ray  fluoroscopic  images  were  acquired  using  an  Acuity 
simulator  (Varian  Medical  Systems,  Palo  Alto,  CA).  The  dimension 
of  the  sensitive  area  of  the  detector  is  397  mm  x  298  mm,  containing 
1024  X  768  pixels.  Two  phantoms  were  used  to  evaluate  the  pre¬ 
sented  KLT  and  PWLS-based  noise  reduction  algorithm.  The  first 
phantom  was  a  commercial  quality  assurance  phantom,  in  which 
several  fine  strips  can  be  used  to  evaluate  the  fine  structure  preserva¬ 
tion  of  the  algorithm.  The  second  phantom  was  an  anthropomorphic 
chest  phantom.  To  simulate  the  respiratory  motion  of  patients,  each 
phantom  was  placed  on  the  top  of  a  platform  capable  of  sinusoidal 
motion  in  the  cranial-caudal  direction.  The  amplitude  of  motion  of 
the  platform  was  3.5  cm,  and  the  period  of  motion  was  set  at  3  s. 
The  selected  motion  parameters  correspond  to  those  of  a  patient 
with  large  respiratory  motion  and  fast  breathing  (22).  In  both  phan¬ 
tom  studies,  the  X-ray  tube  current  was  set  at  10  mA,  and  the  pulse 
length  of  X-ray  was  2  ms  for  low-dose  and  10  ms  for  high-dose  fluo¬ 
roscopic  image.  The  acquisition  frame  rate  was  5  frames  per  second, 
and  all  sequences  consisted  of  20  frames. 


RESULTS 

Quality  assurance  phantom 

We  first  tested  the  presented  algorithm  on  the  quality  as¬ 
surance  (QA)  phantom  placed  on  the  moving  platform. 
The  QA  phantom  contains  several  fine  strips  that  can  be 
used  to  study  the  edge  information  in  the  processed  images. 
Figure  2a  shows  one  frame  of  fluoroscopic  images  obtained 
with  X-ray  tube  current  10  mA  and  duration  of  X-ray  pulse 
length  2  ms.  Figure  2b  shows  the  same  frame  of  image  pro¬ 
cessed  by  the  presented  KL  domain  PWLS  noise  reduction 
algorithm.  It  can  be  observed  that  noise  in  the  low-dose  fluo¬ 
roscopic  image  is  greatly  suppressed,  and  the  image  is  com¬ 
parable  to  that  obtained  with  the  10-ms  protocol  (Fig.  2c). 
One  challenge  for  temporal  filtration  of  fluoroscopic  images 
is  that  tailing  artifacts  may  be  presented  in  the  motion  portion 
of  the  processed  image  (23).  It  can  be  observed  in  Fig.  2b 
that  the  fine  strips  are  well-preserved,  and  no  tailing  artifacts 
have  been  introduced  by  the  temporal  filtering  using  the  KL 
transform.  However,  if  the  isotropic  quadratic  term  was  used 
as  the  penalty  term  in  the  PWLS  objective  function,  the 
edges  in  the  processed  fluoroscopic  images  could  be  blurred. 
Figure  2d  shows  the  image  processed  by  KL  domain  PWLS 
with  isotropic  quadratic  penalty.  It  can  be  observed  that  the 
fine  structures  have  been  severely  blurred.  Profiles  through 
the  fine  strips,  as  shown  in  Fig.  3,  further  illustrate  these  ob¬ 
servations.  Those  results  suggest  that  the  edge-preserving 
penalty  is  desired  in  the  PWLS  objective  function. 
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(d)  (e)  ‘  ‘ 


Fig.  4.  Fluoroscopic  image  of  the  anthropomorphic  chest  phantom,  (a)  Acquired  with  tube  current  10  mA  and  pulse  length 
2  ms;  (b)  image  in  panel  a  processed  by  with  the  two-dimensional  (2D)  mean  filter;  (c)  image  in  panel  a  processed  with  the 
2D  median  filter;  (d)  image  in  panel  a  processed  by  the  Karhunen-Loeve  KL  domain  penalized  weighted  least-squares  noise 
reduction  algorithm;  and  (e)  image  acquired  with  tube  current  10  mA  and  pulse  length  10  ms. 


To  evaluate  the  effectiveness  of  the  presented  algorithm 
quantitatively,  SNR  of  different  regions  of  interest  (ROIs) 
was  calculated.  The  SNR  is  defined  as: 


where  is  the  mean  value  the  signal,  u^,  is  the  mean  value  of 
background  (a  uniform  region  nearby  the  selected  ROI),  and 
is  the  standard  deviation  of  the  signal.  Four  regions  (indi¬ 
cated  by  arrows  in  Fig.  2a)  were  chosen  to  calculate  SNRs, 
and  the  results  are  presented  in  Table  1.  The  improvement  of 
the  SNRs  in  the  image  after  the  KL  domain  PWLS  processing 
varies  according  to  the  locations.  The  range  of  improvement  in 
SNRs  is  between  41%  and  108%,  and  the  average  improve¬ 
ment  is  75%  among  the  four  chosen  ROIs. 

Anthropomorphic  chest  phantom 

Results  of  the  anthropomorphic  chest  are  shown  in  Fig.  4. 
Figure  4a  shows  one  frame  of  fluoroscopic  images  obtained 


with  2-ms  pulse  length  protocol.  Figure  4d  shows  the  image 
processed  by  the  KL  domain  PWLS  noise  reduction  algo¬ 
rithm.  It  can  be  observed  that  the  noise  is  greatly  suppressed, 
whereas  the  fine  structures  are  well  preserved.  Figure  5  pro¬ 
vides  the  difference  image  between  Fig.  4a  and  4d,  in  which 
random  noise  is  dominant  and  no  edges  or  structures  can  be 
observed.  This  indicates  that  good  edge  preservation  can  be 
achieved  by  the  KL  domain  PWLS  noise  reduction  algo¬ 
rithm.  It  can  also  be  observed  that  the  processed  image  ob¬ 
tained  with  the  2-ms  protocol  is  comparable  to  that 
obtained  with  a  10-ms  protocol  (see  Fig.  4e).  We  also  com¬ 
pared  the  presented  algorithm  with  conventional  filters  such 
as  the  mean  and  median  filter.  Figure  4b  shows  the  corre¬ 
sponding  frame  processed  by  the  mean  filter  with  a  3  x  3 
window.  Figure  4c  shows  the  result  from  the  median  filter 
with  a  3  X  3  window.  It  can  be  observed  that  the  image  qual¬ 
ity  of  the  KL  domain  PWLS  processed  image  is  superior  to 
those  images  processed  by  the  mean  and  median  filter  from 
the  bony  stractures  indicated  by  the  rectangular  square  in 
Fig.  4b^d. 
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Fig.  5.  Difference  image  of  the  chest  phantom  between  the  Karhu- 
nen-Loeve  domain  penalized  weighted  least-squares  processed 
image  (Fig.  4d)  and  the  noisy  low-dose  image  (Fig.  4a). 


DISCUSSION 

In  this  work,  we  proposed  a  temporal-spatial  filter  to  reduce 
noise  in  low-dose  X-ray  fluoroscopy.  The  KL  transform  was 
used  to  account  for  correlation  among  neighboring  frames  of 
fluoroscopic  imaging.  The  explicit  estimation  of  motion  fields 
are  avoided  during  the  de-noising  process.  The  KLT  decom¬ 
poses  the  correlated  image  sequences  into  uncorrelated  KL 
components.  Following  KLT,  the  PWLS  criterion  is  used  to  re¬ 
store  each  KL  component,  and  the  smoothing  strength  of  each 
KL  component  is  adaptive  to  its  corresponding  KL  eigenvalue, 
which  reflects  the  SNR  of  the  KL  component.  The  weight  in  the 
PWLS  objective  function  is  determined  by  the  variance  of  each 
measurement,  which  considers  the  signal-dependent  nature  of 
Poisson  noise.  To  preserve  edges  in  the  restored  image,  an  an¬ 
isotropic  quadratic  term  was  used  as  the  penalty  in  the  PWLS 
objective  function.  The  smoothing  strength  of  the  KL  domain 
PWLS  Alter  is  controlled  by  two  parameters:  the  parameter 
6  in  the  weight  coefficient  k'j^,  which  controls  the  strength  of 
local  smoothing  and  the  smoothing  parameter  (S,  which  con¬ 
trols  the  strength  of  global  smoothing. 

In  the  experimental  studies,  the  phantoms  were  placed  on 
a  moving  platform  that  simulates  a  rigid-body  motion.  The 
motion  in  a  real  patient  is  generally  nonrigid  or  deformable. 
The  magnitude  of  deformable  motion  varies  from  voxel  to 
voxel,  whereas  the  motion  in  a  rigid  body  is  uniform  across 


all  voxels.  Thus,  a  deformable  motion  is  generally  more  chal¬ 
lenging  to  handle  with  regard  to  registration  problems  (24). 
However,  the  estimation  of  motion  held  (or  displacement 
field)  between  the  neighboring  images  is  not  required  in  the 
presented  algorithm.  Both  the  similarity  and  the  difference 
are  modeled  through  the  use  of  the  KL  transform.  Because 
of  the  elimination  of  voxel-to-voxel  registration  of  the  neigh¬ 
boring  images,  the  effectiveness  of  the  method  relies  largely 
on  the  behavior  of  the  principle  components  of  the  system.  In 
the  deformable  case,  because  the  motion  occurs  only  on  some 
parts  of  the  image,  the  “average”  motions  are  actually  less 
compared  with  the  rigid-body  motion.  When  the  discrepancy 
between  the  two  images  is  small,  KLT  performance  is  gener¬ 
ally  better  (25).  Therefore,  the  proposed  method  should  work 
as  well  as  it  does  in  the  case  of  rigid-body  motion  investi¬ 
gated  in  this  work.  Indeed,  KLT  has  shown  excellent  perfor¬ 
mance  in  noise  suppression  for  four-dimensional  (4D) 
positron  emission  tomography  (16),  4D  single  photon  emis¬ 
sion  computed  tomography  (25),  and  4D-CT  (26)  in  which 
deformable  motion  exists.  It  is  perhaps  useful  to  mention 
that  a  potential  challenge  when  dealing  with  deformable  mo¬ 
tion  is  the  validation  of  a  fluoroscopic  image  noise  suppres¬ 
sion  algorithm  because  of  the  general  lack  of  the  “ground 
truth.”  This  issue  can  be  partially  resolved  with  the  develop¬ 
ment  of  a  deformable  phantom  (27,  28)  and  the  use  of  inher¬ 
ent  tissue  features  (29). 

In  the  presented  algorithm,  the  KLT  transform  was  used  to 
extract  correlation  information  from  neighboring  images. 
When  neighboring  images  are  not  available  (e.g.,  first  two 
image  slices)  or  the  correlation  between  them  is  weak,  the 
PWLS  criterion  can  be  used  to  suppress  noise  in  the  fluoro¬ 
scopic  image  through  spatial  smoothing.  The  effectiveness 
of  the  PWLS  criterion  in  this  regard  has  been  proven  in  noise 
reduction  of  cone-beam  CT  projection  data  (19,  30)  in  which 
the  data  are  equivalent  to  the  X-ray  fluoroscopic  image.  From 
the  results  in  these  studies  (19,  30),  it  is  reasonable  to  conjec¬ 
ture  that  PWLS  criterion  without  KLT  transform  can  also,  to 
a  certain  extent,  reduce  noise  in  low-dose  X-ray  fluoroscopy 
by  effectively  using  the  spatial  correlation  of  the  voxels 
within  a  fluoroscopic  image. 

In  the  processed  low-dose  fluoroscopic  image  of  two  phan¬ 
toms,  noise  is  greatly  suppressed,  but  fine  stractures  are  well 
preserved.  Comparison  studies  with  traditional  techniques, 
such  as  the  mean  and  median  filters,  show  that  the  proposed 
algorithm  is  advantageous  in  terms  of  structure  preservation. 
The  proposed  noise  reduction  algorithm  can  significantly  im¬ 
prove  the  quality  of  low-dose  X-ray  fluoroscopic  imaging 
and  may  allow  for  image  acquisition  at  significantly  reduced 
doses.  Given  that  frequent  X-ray  imaging  is  increasingly  be¬ 
ing  used  for  therapeutic  guidance,  this  work  should  have  an 
important  impact  on  IGRT  clinical  practice. 
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